diff --git a/Patch_lib/afunction.pro b/Patch_lib/afunction.pro
new file mode 100644
index 0000000000000000000000000000000000000000..19cbd6f2cd440bf256168a9c3c890fada1694d5e
--- /dev/null
+++ b/Patch_lib/afunction.pro
@@ -0,0 +1,90 @@
+function afunction, X, TF
+; afunction - Implements the blurring linear operator via FFT
+; This function applies the blurring operator A to the
+; reconstruction x of the unknown object in the Fourier space, by
+; componentwise multiplication of the two FFTs and then by transforming
+; back the result by an IFFT.
+;
+; SYNOPSIS
+; out = afunction(x, TF)
+;
+; INPUT
+; x (double) - reconstruction
+; TF (dcomplex) - Fourier transform of the blurring operator A
+;
+; OUTPUT
+; out (double) - vectorized blurred version of x
+;
+; ------------------------------------------------------------------------------
+;
+; This software is developed within the research project
+;
+; PRISMA - Optimization methods and software for inverse problems
+; http://www.unife.it/prisma
+;
+; funded by the Italian Ministry for University and Research (MIUR), under
+; the PRIN2008 initiative, grant n. 2008T5KA4L, 2010-2012. This software is
+; part of the package "IRMA - Image Reconstruction in Microscopy and Astronomy"
+; currently under development within the PRISMA project.
+;
+; Version: 1.0
+; Date: July 2011
+
+; Authors:
+; Roberto Cavicchioli, Marco Prato, Luca Zanni
+; Dept. of Pure Appl. Math., Univ. of Modena and Reggio Emilia, Italy
+; roberto.cavicchioli@unimore.it, marco.prato@unimore.it, luca.zanni@unimore.it
+; Mario Bertero, Patrizia Boccacci
+; DISI (Dipartimento di Informatica e Scienze dell'Informazione), University of Genova, Italy
+; bertero@disi.unige.it, boccacci@disi.unige.it
+;
+; Software homepage: http://www.unife.it/prisma/software
+;
+; Copyright (C) 2011 by M. Bertero, P. Boccacci, R. Cavicchioli, M. Prato, L. Zanni.
+; ------------------------------------------------------------------------------
+; COPYRIGHT NOTIFICATION
+;
+; Permission to copy and modify this software and its documentation for
+; internal research use is granted, provided that this notice is retained
+; thereon and on all copies or modifications. The authors and their
+; respective Universities makes no representations as to the suitability
+; and operability of this software for any purpose. It is provided "as is"
+; without express or implied warranty. Use of this software for commercial
+; purposes is expressly prohibited without contacting the authors.
+;
+; This program is free software; you can redistribute it and/or modify it
+; under the terms of the GNU General Public License as published by the
+; Free Software Foundation; either version 3 of the License, or (at your
+; option) any later version.
+;
+; This program is distributed in the hope that it will be useful, but
+; WITHOUT ANY WARRANTY; without even the implied warranty of
+; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+; See the GNU General Public License for more details.
+;
+; You should have received a copy of the GNU General Public License along
+; with this program; if not, either visite http://www.gnu.org/licenses/
+; or write to
+; Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+; ==============================================================================
+
+dimTF = size(TF,/dimensions)
+if (n_elements(dimTF) eq 2) then dimTF = [dimTF,1]
+dimX = size(X,/dimensions)
+out = dblarr(dimTF)
+if (n_elements(dimX) gt 2) then begin
+ for i = 0, dimTF[2]-1 do begin
+ out[*,*,i] = real_part(fft(reform(TF[*,*,i])*fft(reform(X[*,*,i])),/inverse))
+ endfor
+endif else begin
+ x_t = fft(X)
+ for i = 0, dimTF[2]-1 do begin
+ out[*,*,i] = real_part(fft(reform(TF[*,*,i])*x_t,/inverse))
+ endfor
+endelse
+return, out
+
+END
+; ==============================================================================
+; End of afunction.pro file - IRMA package
+; ==============================================================================
diff --git a/Patch_lib/coyote_library/cgscalevector.pro b/Patch_lib/coyote_library/cgscalevector.pro
new file mode 100644
index 0000000000000000000000000000000000000000..75ca92d088793a626a1878df5cf9c796c9b1efc2
--- /dev/null
+++ b/Patch_lib/coyote_library/cgscalevector.pro
@@ -0,0 +1,192 @@
+; docformat = 'rst'
+;
+; NAME:
+; cgScaleVector
+;
+; PURPOSE:
+; This is a utility routine to scale the elements of a vector or an array into a
+; given data range.
+;
+;******************************************************************************************;
+; ;
+; Copyright (c) 1998-2013, by Fanning Software Consulting, Inc. All rights reserved. ;
+; ;
+; Redistribution and use in source and binary forms, with or without ;
+; modification, are permitted provided that the following conditions are met: ;
+; ;
+; * Redistributions of source code must retain the above copyright ;
+; notice, this list of conditions and the following disclaimer. ;
+; * Redistributions in binary form must reproduce the above copyright ;
+; notice, this list of conditions and the following disclaimer in the ;
+; documentation and/or other materials provided with the distribution. ;
+; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ;
+; contributors may be used to endorse or promote products derived from this ;
+; software without specific prior written permission. ;
+; ;
+; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ;
+; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ;
+; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ;
+; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ;
+; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ;
+; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ;
+; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ;
+; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ;
+; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ;
+; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ;
+;******************************************************************************************;
+;
+;+
+; This is a utility routine to scale the elements of a vector or an array into a
+; given data range.
+;
+; :Categories:
+; Utilities
+;
+; :Returns:
+; A vector or array of the same size as the input, scaled into the data range given
+; by `minRange` and `maxRange'. The input vector is confined to the data range set
+; by `MinValue` and `MaxValue` before scaling occurs.
+;
+; :Params:
+; maxRange: in, optional, type=varies, default=1
+; The maximum output value of the scaled vector. Set to 1 by default.
+; minRange: in, optional, type=varies, default=0
+; The minimum output value of the scaled vector. Set to 0 by default.
+; vector: in, required
+; The input vector or array to be scaled.
+;
+; :Keywords:
+; double: in, optional, type=boolean, default=0
+; Set this keyword to perform scaling in double precision. Otherwise, scaling
+; is done in floating point precision.
+; maxvalue: in, optional
+; Set this value to the maximum value of the vector, before scaling (vector < maxvalue).
+; The default value is Max(vector).
+; minvalue: in, optional
+; Set this value to the mimimum value of the vector, before scaling (minvalue < vector).
+; The default value is Min(vector).
+; nan: in, optional, type=boolean, default=0
+; Set this keyword to enable not-a-number checking. NANs in vector will be ignored.
+; preserve_type: in, optional, type=boolean, default=0
+; Set this keyword to preserve the input data type in the output.
+;
+; :Examples:
+; Simple example of scaling a vector::
+;
+; IDL> x = [3, 5, 0, 10]
+; IDL> xscaled = cgScaleVector(x, -50, 50)
+; IDL> Print, xscaled
+; -20.0000 0.000000 -50.0000 50.0000
+
+; Suppose your image has a minimum value of -1.7 and a maximum value = 2.5.
+; You wish to scale this data into the range 0 to 255, but you want to use
+; a diverging color table. Thus, you want to make sure value 0.0 is scaled to 128.
+; You proceed like this::
+;
+; scaledImage = cgScaleVector(image, 0, 255, MINVALUE=-2.5, MAXVALUE=2.5)
+;
+; :Author:
+; FANNING SOFTWARE CONSULTING::
+; David W. Fanning
+; 1645 Sheely Drive
+; Fort Collins, CO 80526 USA
+; Phone: 970-221-0438
+; E-mail: david@idlcoyote.com
+; Coyote's Guide to IDL Programming: http://www.idlcoyote.com
+;
+; :History:
+; Change History::
+; Written by: David W. Fanning, 12 Dec 1998.
+; Added MAXVALUE and MINVALUE keywords. 5 Dec 1999. DWF.
+; Added NAN keyword. 18 Sept 2000. DWF.
+; Removed check that made minRange less than maxRange to allow ranges to be
+; reversed on axes, etc. 28 Dec 2003. DWF.
+; Added PRESERVE_TYPE and DOUBLE keywords. 19 February 2006. DWF.
+; Added FPUFIX to cut down on floating underflow errors. 11 March 2006. DWF.
+; Renamed Scale_Vector to cgScaleVector, 16 May 2013. DWF.
+;
+; :Copyright:
+; Copyright (c) 1998-2013, Fanning Software Consulting, Inc.
+;-
+FUNCTION cgScaleVector, vector, minRange, maxRange, $
+ DOUBLE=double, $
+ MAXVALUE=vectorMax, $
+ MINVALUE=vectorMin, $
+ NAN=nan, $
+ PRESERVE_TYPE=preserve_type
+
+ ; Error handling.
+ Catch, theError
+ IF theError NE 0 THEN BEGIN
+ Catch, /Cancel
+ void = Error_Message()
+ RETURN, vector
+ ENDIF
+
+ ; Check positional parameters.
+ CASE N_Params() OF
+ 0: Message, 'Incorrect number of arguments.'
+ 1: BEGIN
+ IF Keyword_Set(double) THEN BEGIN
+ minRange = 0.0D
+ maxRange = 1.0D
+ ENDIF ELSE BEGIN
+ minRange = 0.0
+ maxRange = 1.0
+ ENDELSE
+ ENDCASE
+ 2: BEGIN
+ IF Keyword_Set(double) THEN maxRange = 1.0D > (minRange + 0.0001D) ELSE $
+ maxRange = 1.0 > (minRange + 0.0001)
+ ENDCASE
+ ELSE:
+ ENDCASE
+
+ ; If input data type is DOUBLE and DOUBLE keyword is not set, then set it.
+ IF Size(FPUFIX(vector), /TNAME) EQ 'DOUBLE' AND N_Elements(double) EQ 0 THEN double = 1
+
+ ; Make sure we are working with at least floating point numbers.
+ IF Keyword_Set(double) THEN minRange = DOUBLE( minRange ) ELSE minRange = FLOAT( minRange )
+ IF Keyword_Set(double) THEN maxRange = DOUBLE( maxRange ) ELSE maxRange = FLOAT( maxRange )
+
+ ; Make sure we have a valid range.
+ IF maxRange EQ minRange THEN Message, 'Range max and min are coincidental'
+
+ ; Check keyword parameters.
+ IF Keyword_Set(double) THEN BEGIN
+ IF N_Elements(vectorMin) EQ 0 THEN vectorMin = Double( Min(FPUFIX(vector), NAN=1) ) $
+ ELSE vectorMin = Double(vectorMin)
+ IF N_Elements(vectorMax) EQ 0 THEN vectorMax = DOUBLE( Max(FPUFIX(vector), NAN=1) ) $
+ ELSE vectorMax = DOUBLE( vectorMax )
+ ENDIF ELSE BEGIN
+ IF N_Elements(vectorMin) EQ 0 THEN vectorMin = FLOAT( Min(FPUFIX(vector), NAN=1) ) $
+ ELSE vectorMin = FLOAT( vectorMin )
+ IF N_Elements(vectorMax) EQ 0 THEN vectorMax = FLOAT( Max(FPUFIX(vector), NAN=Keyword_Set(nan)) ) $
+ ELSE vectorMax = FLOAT( vectorMax )
+ ENDELSE
+
+ ; Trim vector before scaling.
+ index = Where(Finite(vector) EQ 1, count)
+ IF count NE 0 THEN BEGIN
+ IF Keyword_Set(double) THEN trimVector = Double(vector) ELSE trimVector = Float(vector)
+ trimVector[index] = vectorMin > vector[index] < vectorMax
+ ENDIF ELSE BEGIN
+ IF Keyword_Set(double) THEN trimVector = vectorMin > Double(vector) < vectorMax ELSE $
+ trimVector = vectorMin > Float(vector) < vectorMax
+ ENDELSE
+
+ ; Calculate the scaling factors.
+ scaleFactor = [((minRange * vectorMax)-(maxRange * vectorMin)) / $
+ (vectorMax - vectorMin), (maxRange - minRange) / (vectorMax - vectorMin)]
+
+ ; Clear math errors.
+ void = Check_Math()
+
+ ; Return the scaled vector.
+ IF Keyword_Set(preserve_type) THEN BEGIN
+ RETURN, FPUFIX(Convert_To_Type(trimVector * scaleFactor[1] + scaleFactor[0], Size(vector, /TNAME)))
+ ENDIF ELSE BEGIN
+ RETURN, FPUFIX(trimVector * scaleFactor[1] + scaleFactor[0])
+ ENDELSE
+
+END
diff --git a/Patch_lib/coyote_library/error_message.pro b/Patch_lib/coyote_library/error_message.pro
new file mode 100644
index 0000000000000000000000000000000000000000..81dbb37a379369421d23c2041355a1e5b786e6e0
--- /dev/null
+++ b/Patch_lib/coyote_library/error_message.pro
@@ -0,0 +1,280 @@
+;+
+; NAME:
+; ERROR_MESSAGE
+;
+; PURPOSE:
+;
+; The purpose of this function is to have a device-independent
+; error messaging function. The error message is reported
+; to the user by using DIALOG_MESSAGE if widgets are
+; supported and MESSAGE otherwise.
+;
+; In general, the ERROR_MESSAGE function is not called directly.
+; Rather, it is used in a CATCH error handler. Errors are thrown
+; to ERROR_MESSAGE with the MESSAGE command. A typical CATCH error
+; handler is shown below.
+;
+; Catch, theError
+; IF theError NE 0 THEN BEGIN
+; Catch, /Cancel
+; void = Error_Message()
+; RETURN
+; ENDIF
+;
+; Error messages would get into the ERROR_MESSAGE function by
+; throwing an error with the MESSAGE command, like this:
+;
+; IF test NE 1 THEN Message, 'The test failed.'
+;
+; AUTHOR:
+;
+; FANNING SOFTWARE CONSULTING
+; David Fanning, Ph.D.
+; 1645 Sheely Drive
+; Fort Collins, CO 80526 USA
+; Phone: 970-221-0438
+; E-mail: david@idlcoyote.com
+; Coyote's Guide to IDL Programming: http://www.idlcoyote.com/
+;
+; CATEGORY:
+;
+; Utility.
+;
+; CALLING SEQUENCE:
+;
+; ok = Error_Message(the_Error_Message)
+;
+; INPUTS:
+;
+; the_Error_Message: This is a string argument containing the error
+; message you want reported. If undefined, this variable is set
+; to the string in the !Error_State.Msg system variable.
+;
+; KEYWORDS:
+;
+; ERROR: Set this keyword to cause Dialog_Message to use the ERROR
+; reporting dialog. Note that a bug in IDL causes the ERROR dialog
+; to be used whether this keyword is set to 0 or 1!
+;
+; INFORMATIONAL: Set this keyword to cause Dialog_Message to use the
+; INFORMATION dialog instead of the WARNING dialog. Note that a bug
+; in IDL causes the ERROR dialog to be used if this keyword is set to 0!
+;
+; NONAME: Normally, the name of the routine in which the error occurs is
+; added to the error message. Setting this keyword will suppress this
+; behavior.
+;
+; TITLE: Set this keyword to the title of the DIALOG_MESSAGE window. By
+; default the keyword is set to 'System Error' unless !ERROR_STATE.NAME
+; equals "IDL_M_USER_ERR", in which case it is set to "Trapped Error'.
+;
+; TRACEBACK: Setting this keyword results in an error traceback
+; being printed to standard output with the PRINT command. Set to
+; 1 (ON) by default. Use TRACEBACK=0 to turn this functionality off.
+;
+; QUIET: Set this keyword to suppress the DIALOG_MESSAGE pop-up dialog.
+;
+; OUTPUTS:
+;
+; Currently the only output from the function is the string "OK".
+;
+; RESTRICTIONS:
+;
+; The WARNING Dialog_Message dialog is used by default.
+;
+; EXAMPLE:
+;
+; To handle an undefined variable error:
+;
+; IF N_Elements(variable) EQ 0 THEN $
+; ok = Error_Message('Variable is undefined')
+;
+; MODIFICATION HISTORY:
+;
+; Written by: David W. Fanning, 27 April 1999.
+; Added the calling routine's name in the message and NoName keyword. 31 Jan 2000. DWF.
+; Added _Extra keyword. 10 February 2000. DWF.
+; Forgot to add _Extra everywhere. Fixed for MAIN errors. 8 AUG 2000. DWF.
+; Adding call routine's name to Traceback Report. 8 AUG 2000. DWF.
+; Added ERROR, INFORMATIONAL, and TITLE keywords. 19 SEP 2002. DWF.
+; Removed the requirement that you use the NONAME keyword with the MESSAGE
+; command when generating user-trapped errors. 19 SEP 2002. DWF.
+; Added distinctions between trapped errors (errors generated with the
+; MESSAGE command) and IDL system errors. Note that if you call ERROR_MESSAGE
+; directly, then the state of the !ERROR_STATE.NAME variable is set
+; to the *last* error generated. It is better to access ERROR_MESSAGE
+; indirectly in a Catch error handler from the MESSAGE command. 19 SEP 2002. DWF.
+; Change on 19 SEP 2002 to eliminate NONAME requirement did not apply to object methods.
+; Fixed program to also handle messages from object methods. 30 JULY 2003. DWF.
+; Removed obsolete STR_SEP and replaced with STRSPLIT. 27 Oct 2004. DWF.
+; Made a traceback the default case without setting TRACEBACK keyword. 19 Nov 2004. DWF.
+; Added check for window connection specifically for CRON jobs. 6 May 2008. DWF.
+; Added QUIET keyword. 18 October 2008. DWF.
+; The traceback information was bypassed when in the PostScript device. Not what I
+; had in mind. Fixed. 6 July 2009. DWF.
+; The QUIET keyword was clearing traceback information. Fixed with help from Phillip Bitzer. 2 Oct 2012. DWF.
+;-
+;******************************************************************************************;
+; Copyright (c) 2008, by Fanning Software Consulting, Inc. ;
+; All rights reserved. ;
+; ;
+; Redistribution and use in source and binary forms, with or without ;
+; modification, are permitted provided that the following conditions are met: ;
+; ;
+; * Redistributions of source code must retain the above copyright ;
+; notice, this list of conditions and the following disclaimer. ;
+; * Redistributions in binary form must reproduce the above copyright ;
+; notice, this list of conditions and the following disclaimer in the ;
+; documentation and/or other materials provided with the distribution. ;
+; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ;
+; contributors may be used to endorse or promote products derived from this ;
+; software without specific prior written permission. ;
+; ;
+; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ;
+; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ;
+; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ;
+; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ;
+; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ;
+; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ;
+; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ;
+; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ;
+; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ;
+; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ;
+;******************************************************************************************;
+FUNCTION ERROR_MESSAGE, theMessage, Error=error, Informational=information, $
+ Traceback=traceback, NoName=noname, Title=title, Quiet=quiet, _Extra=extra
+
+ On_Error, 2
+
+ ; Check for presence and type of message.
+
+ IF N_Elements(theMessage) EQ 0 THEN theMessage = !Error_State.Msg
+ s = Size(theMessage)
+ messageType = s[s[0]+1]
+ IF messageType NE 7 THEN BEGIN
+ Message, "The message parameter must be a string.", _Extra=extra
+ ENDIF
+ IF N_Elements(traceback) EQ 0 THEN traceback = 1
+
+ ; Get the call stack and the calling routine's name.
+ Help, Calls=callStack
+ callingRoutine = (StrSplit(StrCompress(callStack[1])," ", /Extract))[0]
+
+ ; Are widgets supported?
+ IF !D.Name EQ 'PS' OR !D.Name EQ 'Z' THEN BEGIN
+ widgetsSupported = 1
+ ENDIF ELSE BEGIN
+ widgetsSupported = ((!D.Flags AND 65536L) NE 0)
+ ENDELSE
+
+ ; Is the QUIET keyword set? Then no dialogs.
+ IF Keyword_Set(quiet) THEN widgetsSupported = 0
+
+ ; It is not enough to know if widgets are supported. In CRON jobs, widgets are
+ ; supported, but there is no X connection and pop-up dialogs are not allowed.
+ ; Here is a quick test to see if we can connect to a windowing system. If not,
+ ; then we are going to assume widgets are not supported.
+ Catch, theError
+ IF theError NE 0 THEN BEGIN
+ Catch, /CANCEL
+ widgetsSupported = 0
+ GOTO, testWidgetSupport
+ ENDIF
+ theWindow = !D.Window
+ IF (!D.Flags AND 256) NE 0 THEN Window, /FREE, XSIZE=5, YSIZE=5, /PIXMAP
+ Catch, /CANCEL
+
+ testWidgetSupport: ; Come here if you choke on creating a window.
+ IF !D.Window NE theWindow THEN BEGIN
+ WDelete, !D.Window
+ IF theWindow GE 0 THEN WSet, theWindow
+ ENDIF
+
+ IF widgetsSupported THEN BEGIN
+
+ ; If this is an error produced with the MESSAGE command, it is a trapped
+ ; error and will have the name "IDL_M_USER_ERR".
+ IF !ERROR_STATE.NAME EQ "IDL_M_USER_ERR" THEN BEGIN
+
+ IF N_Elements(title) EQ 0 THEN title = 'Trapped Error'
+
+ ; If the message has the name of the calling routine in it,
+ ; it should be stripped out. Can you find a colon in the string?
+
+ ; Is the calling routine an object method? If so, special processing
+ ; is required. Object methods will have two colons together.
+ doublecolon = StrPos(theMessage, "::")
+ IF doublecolon NE -1 THEN BEGIN
+
+ prefix = StrMid(theMessage, 0, doublecolon+2)
+ submessage = StrMid(theMessage, doublecolon+2)
+ colon = StrPos(submessage, ":")
+ IF colon NE -1 THEN BEGIN
+
+ ; Extract the text up to the colon. Is this the same as
+ ; the callingRoutine? If so, strip it.
+ IF StrMid(theMessage, 0, colon+StrLen(prefix)) EQ callingRoutine THEN $
+ theMessage = StrMid(theMessage, colon+1+StrLen(prefix))
+ ENDIF
+ ENDIF ELSE BEGIN
+
+ colon = StrPos(theMessage, ":")
+ IF colon NE -1 THEN BEGIN
+
+ ; Extract the text up to the colon. Is this the same as
+ ; the callingRoutine? If so, strip it.
+ IF StrMid(theMessage, 0, colon) EQ callingRoutine THEN $
+ theMessage = StrMid(theMessage, colon+1)
+ ENDIF
+
+ ENDELSE
+
+
+ ; Add the calling routine's name, unless NONAME is set.
+ IF Keyword_Set(noname) THEN BEGIN
+ answer = Dialog_Message(theMessage, Title=title, _Extra=extra, $
+ Error=error, Information=information)
+ ENDIF ELSE BEGIN
+ answer = Dialog_Message(StrUpCase(callingRoutine) + ": " + $
+ theMessage, Title=title, _Extra=extra, $
+ Error=error, Information=information)
+ ENDELSE
+
+ ENDIF ELSE BEGIN
+
+ ; Otherwise, this is an IDL system error.
+ IF N_Elements(title) EQ 0 THEN title = 'System Error'
+
+ IF StrUpCase(callingRoutine) EQ "$MAIN$" THEN $
+ answer = Dialog_Message(theMessage, _Extra=extra, Title=title, $
+ Error=error, Information=information) ELSE $
+ IF Keyword_Set(noname) THEN BEGIN
+ answer = Dialog_Message(theMessage, _Extra=extra, Title=title, $
+ Error=error, Information=information)
+ ENDIF ELSE BEGIN
+ answer = Dialog_Message(StrUpCase(callingRoutine) + "--> " + $
+ theMessage, _Extra=extra, Title=title, $
+ Error=error, Information=information)
+ ENDELSE
+ ENDELSE
+ ENDIF ELSE BEGIN
+ Help, /Last_Message, Output=traceback_msg ; Required because following MESSAGE call clears traceback info.
+ Message, theMessage, /Continue, /NoPrint, /NoName, /NoPrefix, _Extra=extra
+ IF Keyword_Set(noname) THEN $
+ Print, theMessage ELSE $
+ Print, '%' + callingRoutine + ': ' + theMessage
+ answer = 'OK'
+ ENDELSE
+
+ ; Provide traceback information if requested and this is NOT an informational message.
+ IF Keyword_Set(traceback) AND ~Keyword_Set(informational)THEN BEGIN
+ IF N_Elements(traceback_msg) NE 0 THEN traceback = traceback_msg ELSE Help, /Last_Message, Output=traceback
+ Print,''
+ Print, 'Traceback Report from ' + StrUpCase(callingRoutine) + ':'
+ Print, ''
+ FOR j=0,N_Elements(traceback)-1 DO Print, " " + traceback[j]
+ ENDIF
+
+ RETURN, answer
+END ; ----------------------------------------------------------------------------
+
diff --git a/Patch_lib/coyote_library/fpufix.pro b/Patch_lib/coyote_library/fpufix.pro
new file mode 100644
index 0000000000000000000000000000000000000000..d958b61c9a955ee2342672d705edea347b71c47a
--- /dev/null
+++ b/Patch_lib/coyote_library/fpufix.pro
@@ -0,0 +1,136 @@
+;+
+; NAME:
+; FPUFIX
+;
+; PURPOSE:
+;
+; This is a utility routine to examine a variable and fix problems
+; that will create floating point underflow errors.
+;
+; AUTHOR:
+;
+; FANNING SOFTWARE CONSULTING
+; David Fanning, Ph.D.
+; 1645 Sheely Drive
+; Fort Collins, CO 80526 USA
+; Phone: 970-221-0438
+; E-mail: david@idlcoyote.com
+; Coyote's Guide to IDL Programming: http://www.idlcoyote.com
+;
+; CATEGORY:
+;
+; Utilities
+;
+; CALLING SEQUENCE:
+;
+; fixedData = FPUFIX(data)
+;
+; ARGUMENTS:
+;
+; data : A numerical variable to be checked for values that will cause
+; floating point underflow errors. Suspect values are set to 0.
+;
+; KEYWORDS:
+;
+; None.
+
+; RETURN VALUE:
+;
+; fixedData: The output is the same as the input, except that any values that
+; will cause subsequent floating point underflow errors are set to 0.
+;
+; COMMON BLOCKS:
+; None.
+;
+; EXAMPLES:
+;
+; data = FPTFIX(data)
+;
+; RESTRICTIONS:
+;
+; None.
+;
+; MODIFICATION HISTORY:
+;
+; Written by David W. Fanning, from Mati Meron's example FPU_FIX. Mati's
+; program is more robust that this (ftp://cars3.uchicago.edu/midl/),
+; but this serves my needs and doesn't require other programs from
+; Mati's library. 24 February 2006.
+;-
+;******************************************************************************************;
+; Copyright (c) 2008, by Fanning Software Consulting, Inc. ;
+; All rights reserved. ;
+; ;
+; Redistribution and use in source and binary forms, with or without ;
+; modification, are permitted provided that the following conditions are met: ;
+; ;
+; * Redistributions of source code must retain the above copyright ;
+; notice, this list of conditions and the following disclaimer. ;
+; * Redistributions in binary form must reproduce the above copyright ;
+; notice, this list of conditions and the following disclaimer in the ;
+; documentation and/or other materials provided with the distribution. ;
+; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ;
+; contributors may be used to endorse or promote products derived from this ;
+; software without specific prior written permission. ;
+; ;
+; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ;
+; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ;
+; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ;
+; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ;
+; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ;
+; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ;
+; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ;
+; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ;
+; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ;
+; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ;
+;******************************************************************************************;
+FUNCTION FPUFIX, data
+
+ ; Return to caller on error after setting !Except.
+ On_Error, 2
+ Catch, theError
+ IF theError NE 0 THEN BEGIN
+ Catch, /Cancel
+ Message, !Error_State.Msg
+ IF N_Elements(except) THEN !EXCEPT = except
+ RETURN, data
+ ENDIF
+
+ ; Only want to deal with numerical data types.
+ ; Return all other kinds.
+ dataType = Size(data, /Type)
+ nogoodtypes = [0,7,8,10,11]
+ void = Where(nogoodTypes EQ dataType, count)
+ IF count GT 0 THEN RETURN, data
+
+ ; Floating underflow error we are trying to fix.
+ fpu_error = 32
+
+ ; Save current !EXCEPT. Don't report exceptions here.
+ except = !EXCEPT
+ !EXCEPT = 1
+
+ ; Clear math error status.
+ void = Check_Math()
+
+ ; Do something with the data that will cause floating underflow errors.
+ void = Min(data, /NAN)
+
+ ; Check the math error status now.
+ check = Check_Math()
+
+ ; If this is a floating underflow error, then fix it.
+ IF check EQ fpu_error THEN BEGIN
+ info = MaChar(DOUBLE=(dataType EQ 5 OR dataType EQ 9))
+ indices = Where(Abs(data) LT info.xmin, count)
+ IF count GT 0 THEN data[indices] = 0
+ ENDIF
+
+ ; Clean up.
+ !EXCEPT = except
+ void = Check_Math()
+
+ ; Return the repaired data.
+ RETURN, data
+
+END
diff --git a/Patch_lib/coyote_library/gmascl.pro b/Patch_lib/coyote_library/gmascl.pro
new file mode 100644
index 0000000000000000000000000000000000000000..70bc36d20aeac32427c5f71c20097521de046c7a
--- /dev/null
+++ b/Patch_lib/coyote_library/gmascl.pro
@@ -0,0 +1,166 @@
+;+
+; NAME:
+; GMASCL
+;
+; PURPOSE:
+;
+; This is a utility routine to perform basic gray-level pixel
+; transformations of images. I think of it as BYTSCL on steroids.
+; It is similar to IMADJUST in _Digital Image Processing with MATLAB_
+; by Gonzales, Wood, and Eddins. It performs a log scaling of the
+; image array.
+;
+; AUTHOR:
+;
+; FANNING SOFTWARE CONSULTING
+; David Fanning, Ph.D.
+; 1645 Sheely Drive
+; Fort Collins, CO 80526 USA
+; Phone: 970-221-0438
+; E-mail: david@idlcoyote.com
+; Coyote's Guide to IDL Programming: http://www.idlcoyote.com
+;
+; CATEGORY:
+;
+; Utilities
+;
+; CALLING SEQUENCE:
+;
+; scaledImage = GMASCL(image)
+;
+; ARGUMENTS:
+;
+; image: The image to be scaled. Written for 2D images, but arrays
+; of any size are treated alike.
+;
+; KEYWORDS:
+;
+; GAMMA: The exponent in a power-law transformation (image^gamma). A gamma
+; value of 1 results in a linear distribution of values between
+; OMIN and OMAX. Gamma values less than 1 map darker image values
+; into a wider range of output values, and Gamma values
+; greater than 1 maps lighter image values into a wider
+; range of output values. The gamma value is constrained to
+; be greater than 1.0e-6.
+;
+; MAX: Any value in the input image greater than this value is
+; set to this value before scaling.
+;
+; MIN: Any value in the input image less than this value is
+; set to this value before scaling.
+;
+; NEGATIVE: If set, the "negative" of the result is returned.
+;
+; OMAX: The output image is scaled between OMIN and OMAX. The
+; default value is 255.
+;
+; OMIN: The output image is scaled between OMIN and OMAX. The
+; default value is 0.
+; RETURN VALUE:
+;
+; scaledImage: The output, scaled into the range OMIN to OMAX. A byte array.
+;
+; COMMON BLOCKS:
+; None.
+;
+; EXAMPLES:
+;
+; LoadCT, 0 ; Gray-scale colors.
+; image = cgDemoData(11) ; Load image.
+; TV, GmaScl(image, Min=30, Max=100) ; Similar to BytScl.
+; TV, GmaScl(image, /Negative) ; Produce negative image.
+; power = Shift(ALog(Abs(FFT(image,-1))), 124, 124) ; Create power spectrum.
+; TV, GmaScl(power, Gamma=2.5) ; View power specturm with gamma correction.
+; TV, GmaScl(power, Gamma=2.5, /Negative) ; Reverse power spectrum.
+;
+; RESTRICTIONS:
+;
+; Requires cgScaleVector from the Coyote Library:
+;
+; http://www.idlcoyote.com/programs/cgScaleVector.pro
+;
+; MODIFICATION HISTORY:
+;
+; Written by: David W. Fanning, 17 February 2006.
+; Fixed a problem with output scaling. 1 July 2009. DWF (with input from Bo Milvang-Jensen).
+; Now setting NAN keyword on all MIN and MAX functions. 2 Dec 2011. DWF.
+;-
+;******************************************************************************************;
+; Copyright (c) 2008, by Fanning Software Consulting, Inc. ;
+; All rights reserved. ;
+; ;
+; Redistribution and use in source and binary forms, with or without ;
+; modification, are permitted provided that the following conditions are met: ;
+; ;
+; * Redistributions of source code must retain the above copyright ;
+; notice, this list of conditions and the following disclaimer. ;
+; * Redistributions in binary form must reproduce the above copyright ;
+; notice, this list of conditions and the following disclaimer in the ;
+; documentation and/or other materials provided with the distribution. ;
+; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ;
+; contributors may be used to endorse or promote products derived from this ;
+; software without specific prior written permission. ;
+; ;
+; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ;
+; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ;
+; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ;
+; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ;
+; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ;
+; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ;
+; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ;
+; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ;
+; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ;
+; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ;
+;******************************************************************************************;
+FUNCTION GmaScl, image, $
+ GAMMA=gamma, $
+ MAX=imageMax, $
+ MIN=imageMin, $
+ NEGATIVE=negative, $
+ OMAX=maxOut, $
+ OMIN=minOut
+
+ ; Return to caller on error.
+ ;On_Error, 2
+ Catch, theError
+ IF theError NE 0 THEN BEGIN
+ Catch, /Cancel
+ void = Error_Message()
+ RETURN, vector
+ ENDIF
+
+ ; Check arguments.
+ IF N_Elements(image) EQ 0 THEN Message, 'Must pass IMAGE argument.'
+
+ ; Check for underflow of values near 0. Yuck!
+ curExcept = !Except
+ !Except = 0
+ i = Where(image GT -1e-35 AND image LT 1e-35, count)
+ IF count GT 0 THEN image[i] = 0.0
+ void = Check_Math()
+ !Except = curExcept
+
+ output = Double(image)
+
+ ; Check keywords.
+ IF N_Elements(imageMax) EQ 0 THEN imageMax = Max(output, /NAN)
+ IF N_Elements(imageMin) EQ 0 THEN imageMin = Min(output, /NAN)
+ IF N_Elements(maxOut) EQ 0 THEN maxOut = 255B ELSE maxout = 0 > Byte(maxOut) < 255
+ IF N_Elements(minOut) EQ 0 THEN minOut = 0B ELSE minOut = 0 > Byte(minOut) < 255
+ IF minOut GE maxout THEN Message, 'OMIN must be less than OMAX.'
+
+ ; Gamma must be greater than 0.
+ IF N_Elements(gamma) EQ 0 THEN gamma = 1.0D ELSE gamma = 1.0e-6 > Double(gamma)
+
+ ; Perform initial scaling of the image.
+ output = cgScaleVector(Temporary(output), 0.0D, 1.0D, MinValue=imageMin, MaxValue=imageMax, /NAN, Double=1)
+
+ ; For gamma, we need positive values. Make sure we have them
+ IF Min(output, /NAN) LT 0.0D THEN output = output + Abs(Min(output, /NAN))
+ output = cgScaleVector(output^gamma, minOut, maxOut, MinValue=0.0D, MaxValue=1.0D, /NAN, Double=1)
+
+ ; Does the user want the negative result?
+ IF Keyword_Set(negative) THEN RETURN, BYTE(0B > (maxout - Round(output) + minOut) < 255B) $
+ ELSE RETURN, BYTE(0B > Round(output) < 255B)
+
+END
diff --git a/Patch_lib/coyote_library/logscl.pro b/Patch_lib/coyote_library/logscl.pro
new file mode 100644
index 0000000000000000000000000000000000000000..fbd9aab9b619dcbc561b63c27ef1d6a1fd508802
--- /dev/null
+++ b/Patch_lib/coyote_library/logscl.pro
@@ -0,0 +1,176 @@
+;+
+; NAME:
+; LOGSCL
+;
+; PURPOSE:
+;
+; This is a utility routine to perform a log intensity transformation
+; on an image. For exponent values greater than 1.0, the upper and
+; lower values of the image are compressed and centered on the mean.
+; Larger exponent values provide steeper compression. For exponent values
+; less than 1.0, the compression is similar to gamma compression. (See
+; IMGSCL.) See pages 68-70 in _Digital Image Processing with MATLAB_
+; by Gonzales, Wood, and Eddins. The function is used to improve contrast
+; in images.
+;
+; AUTHOR:
+;
+; FANNING SOFTWARE CONSULTING
+; David Fanning, Ph.D.
+; 1645 Sheely Drive
+; Fort Collins, CO 80526 USA
+; Phone: 970-221-0438
+; E-mail: david@idlcoyote.com
+; Coyote's Guide to IDL Programming: http://www.idlcoyote.com
+;
+; CATEGORY:
+;
+; Utilities
+;
+; CALLING SEQUENCE:
+;
+; outputImage = LOGSCL(image)
+;
+; ARGUMENTS:
+;
+; image: The image to be scaled. Written for 2D images, but arrays
+; of any size are treated alike.
+;
+; KEYWORDS:
+;
+; EXPONENT: The exponent in a log transformation. By default, 4.0.
+;
+; MEAN: Values on either side of the mean will be compressed by the log.
+; The value is a normalized value between 0.0 and 1.0. By default, 0.5.
+;
+; NEGATIVE: If set, the "negative" of the result is returned.
+;
+; MAX: Any value in the input image greater than this value is
+; set to this value before scaling.
+;
+; MIN: Any value in the input image less than this value is
+; set to this value before scaling.
+;
+; OMAX: The output image is scaled between OMIN and OMAX. The
+; default value is 255.
+;
+; OMIN: The output image is scaled between OMIN and OMAX. The
+; default value is 0.
+; RETURN VALUE:
+;
+; outputImage: The output, scaled into the range OMIN to OMAX. A byte array.
+;
+; COMMON BLOCKS:
+; None.
+;
+; EXAMPLES:
+;
+; cgLoadCT, 0 ; Gray-scale colors.
+; image = cgDemoData(22) ; Load image.
+; cgImage, image ; No contrast.
+; cgImage, LogScl(image) ; Improved contrast.
+; cgImage, LogScl(image, Exponent=10, Mean=0.65) ; Even more contrast.
+; cgImage, LogScl(image, /Negative, Exponent=5) ; A negative image.
+;
+; RESTRICTIONS:
+;
+; Requires cgScaleVector from the Coyote Library:
+;
+; http://www.idlcoyote.com/programs/cgScaleVector.pro
+;
+; MODIFICATION HISTORY:
+;
+; Written by: David W. Fanning, 20 February 2006.
+; Fixed a problem with output scaling. 1 July 2009. DWF (with input from Bo Milvang-Jensen).
+;-
+;******************************************************************************************;
+; Copyright (c) 2008, by Fanning Software Consulting, Inc. ;
+; All rights reserved. ;
+; ;
+; Redistribution and use in source and binary forms, with or without ;
+; modification, are permitted provided that the following conditions are met: ;
+; ;
+; * Redistributions of source code must retain the above copyright ;
+; notice, this list of conditions and the following disclaimer. ;
+; * Redistributions in binary form must reproduce the above copyright ;
+; notice, this list of conditions and the following disclaimer in the ;
+; documentation and/or other materials provided with the distribution. ;
+; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ;
+; contributors may be used to endorse or promote products derived from this ;
+; software without specific prior written permission. ;
+; ;
+; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ;
+; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ;
+; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ;
+; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ;
+; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ;
+; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ;
+; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ;
+; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ;
+; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ;
+; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ;
+;******************************************************************************************;
+FUNCTION LogScl, image, $
+ EXPONENT=exponent, $
+ MEAN=mean, $
+ NEGATIVE=negative, $
+ MAX=maxValue, $
+ MIN=minValue, $
+ OMAX=maxOut, $
+ OMIN=minOut
+
+ ; Return to caller on error.
+ ;On_Error, 2
+ Catch, theError
+ IF theError NE 0 THEN BEGIN
+ Catch, /Cancel
+ void = Error_Message()
+ RETURN, vector
+ ENDIF
+
+ ; Check arguments.
+ IF N_Elements(image) EQ 0 THEN Message, 'Must pass IMAGE argument.'
+
+ ; Check for underflow of values near 0. Yuck!
+ curExcept = !Except
+ !Except = 0
+ i = Where(image GT -1e-35 AND image LT 1e-35, count)
+ IF count GT 0 THEN image[i] = 0.0
+ void = Check_Math()
+ !Except = curExcept
+
+ output = Double(image)
+
+ ; Check keywords.
+ IF N_Elements(exponent) EQ 0 THEN exponent = 4.0
+ IF N_Elements(mean) EQ 0 THEN mean = 0.5 ELSE mean = 0.0 > mean < 1.0
+ IF N_Elements(maxOut) EQ 0 THEN maxOut = 255B ELSE maxout = 0 > Byte(maxOut) < 255
+ IF N_Elements(minOut) EQ 0 THEN minOut = 0B ELSE minOut = 0 > Byte(minOut) < 255
+ IF minOut GE maxout THEN Message, 'OMIN must be less than OMAX.'
+
+ ; Exponent must be greater than 0.
+ exponent = 1.0e-6 > exponent
+
+ ; Perform initial scaling of the image into 0 to 1.
+ output = cgScaleVector(Temporary(output), 0.0D, 1.0D, MaxValue=maxValue, $
+ MinValue=minValue, /NAN, Double=1)
+
+ ; Too damn many floating underflow warnings, no matter WHAT I do! :-(
+ thisExcept = !Except
+ !Except = 0
+
+ ; Log scaling.
+ output = 1.0D / ((1.0D + (mean / (Temporary(output) > 1e-16))^exponent) > (1e-16))
+
+ ; Scale to image coordinates.
+ output = cgScaleVector(Temporary(output), minOut, maxOut, MinValue=0.0D, MaxValue=1.0D, /NAN, Double=1)
+
+ ; Clear math errors.
+ void = Check_Math()
+ !Except = thisExcept
+
+ ; Does the user want the negative result?
+ IF Keyword_Set(negative) THEN RETURN, BYTE(maxout - Round(output) + minOut) $
+ ELSE RETURN, BYTE(Round(output))
+
+END
diff --git a/Patch_lib/coyote_library/scale_vector.pro b/Patch_lib/coyote_library/scale_vector.pro
new file mode 100644
index 0000000000000000000000000000000000000000..2776b3a6d92a0b2ee2c7af36bf4394ec65b2af6b
--- /dev/null
+++ b/Patch_lib/coyote_library/scale_vector.pro
@@ -0,0 +1,192 @@
+; docformat = 'rst'
+;
+; NAME:
+; Scale_Vector
+;
+; PURPOSE:
+; This is a utility routine to scale the elements of a vector or an array into a
+; given data range.
+;
+;******************************************************************************************;
+; ;
+; Copyright (c) 1998-2013, by Fanning Software Consulting, Inc. All rights reserved. ;
+; ;
+; Redistribution and use in source and binary forms, with or without ;
+; modification, are permitted provided that the following conditions are met: ;
+; ;
+; * Redistributions of source code must retain the above copyright ;
+; notice, this list of conditions and the following disclaimer. ;
+; * Redistributions in binary form must reproduce the above copyright ;
+; notice, this list of conditions and the following disclaimer in the ;
+; documentation and/or other materials provided with the distribution. ;
+; * Neither the name of Fanning Software Consulting, Inc. nor the names of its ;
+; contributors may be used to endorse or promote products derived from this ;
+; software without specific prior written permission. ;
+; ;
+; THIS SOFTWARE IS PROVIDED BY FANNING SOFTWARE CONSULTING, INC. ''AS IS'' AND ANY ;
+; EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES ;
+; OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT ;
+; SHALL FANNING SOFTWARE CONSULTING, INC. BE LIABLE FOR ANY DIRECT, INDIRECT, ;
+; INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED ;
+; TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; ;
+; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ;
+; ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT ;
+; (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS ;
+; SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. ;
+;******************************************************************************************;
+;
+;+
+; This is a utility routine to scale the elements of a vector or an array into a
+; given data range. The routine is now depreciated in favor of cgScaleVector.
+;
+; :Categories:
+; Utilities
+;
+; :Returns:
+; A vector or array of the same size as the input, scaled into the data range given
+; by `minRange` and `maxRange'. The input vector is confined to the data range set
+; by `MinValue` and `MaxValue` before scaling occurs.
+;
+; :Params:
+; maxRange: in, optional, type=varies, default=1
+; The maximum output value of the scaled vector. Set to 1 by default.
+; minRange: in, optional, type=varies, default=0
+; The minimum output value of the scaled vector. Set to 0 by default.
+; vector: in, required
+; The input vector or array to be scaled.
+;
+; :Keywords:
+; double: in, optional, type=boolean, default=0
+; Set this keyword to perform scaling in double precision. Otherwise, scaling
+; is done in floating point precision.
+; maxvalue: in, optional
+; Set this value to the maximum value of the vector, before scaling (vector < maxvalue).
+; The default value is Max(vector).
+; minvalue: in, optional
+; Set this value to the mimimum value of the vector, before scaling (minvalue < vector).
+; The default value is Min(vector).
+; nan: in, optional, type=boolean, default=0
+; Set this keyword to enable not-a-number checking. NANs in vector will be ignored.
+; preserve_type: in, optional, type=boolean, default=0
+; Set this keyword to preserve the input data type in the output.
+;
+; :Examples:
+; Simple example of scaling a vector::
+;
+; IDL> x = [3, 5, 0, 10]
+; IDL> xscaled = Scale_Vector(x, -50, 50)
+; IDL> Print, xscaled
+; -20.0000 0.000000 -50.0000 50.0000
+
+; Suppose your image has a minimum value of -1.7 and a maximum value = 2.5.
+; You wish to scale this data into the range 0 to 255, but you want to use
+; a diverging color table. Thus, you want to make sure value 0.0 is scaled to 128.
+; You proceed like this::
+;
+; scaledImage = Scale_Vector(image, 0, 255, MINVALUE=-2.5, MAXVALUE=2.5)
+;
+; :Author:
+; FANNING SOFTWARE CONSULTING::
+; David W. Fanning
+; 1645 Sheely Drive
+; Fort Collins, CO 80526 USA
+; Phone: 970-221-0438
+; E-mail: david@idlcoyote.com
+; Coyote's Guide to IDL Programming: http://www.idlcoyote.com
+;
+; :History:
+; Change History::
+; Written by: David W. Fanning, 12 Dec 1998.
+; Added MAXVALUE and MINVALUE keywords. 5 Dec 1999. DWF.
+; Added NAN keyword. 18 Sept 2000. DWF.
+; Removed check that made minRange less than maxRange to allow ranges to be
+; reversed on axes, etc. 28 Dec 2003. DWF.
+; Added PRESERVE_TYPE and DOUBLE keywords. 19 February 2006. DWF.
+; Added FPUFIX to cut down on floating underflow errors. 11 March 2006. DWF.
+; This routine is being depreciated in favor of cgScaleVector. 16 May 2013. DWF.
+;
+; :Copyright:
+; Copyright (c) 1998-2013, Fanning Software Consulting, Inc.
+;-
+FUNCTION Scale_Vector, vector, minRange, maxRange, $
+ DOUBLE=double, $
+ MAXVALUE=vectorMax, $
+ MINVALUE=vectorMin, $
+ NAN=nan, $
+ PRESERVE_TYPE=preserve_type
+
+ ; Error handling.
+ Catch, theError
+ IF theError NE 0 THEN BEGIN
+ Catch, /Cancel
+ void = Error_Message()
+ RETURN, vector
+ ENDIF
+
+ ; Check positional parameters.
+ CASE N_Params() OF
+ 0: Message, 'Incorrect number of arguments.'
+ 1: BEGIN
+ IF Keyword_Set(double) THEN BEGIN
+ minRange = 0.0D
+ maxRange = 1.0D
+ ENDIF ELSE BEGIN
+ minRange = 0.0
+ maxRange = 1.0
+ ENDELSE
+ ENDCASE
+ 2: BEGIN
+ IF Keyword_Set(double) THEN maxRange = 1.0D > (minRange + 0.0001D) ELSE $
+ maxRange = 1.0 > (minRange + 0.0001)
+ ENDCASE
+ ELSE:
+ ENDCASE
+
+ ; If input data type is DOUBLE and DOUBLE keyword is not set, then set it.
+ IF Size(FPUFIX(vector), /TNAME) EQ 'DOUBLE' AND N_Elements(double) EQ 0 THEN double = 1
+
+ ; Make sure we are working with at least floating point numbers.
+ IF Keyword_Set(double) THEN minRange = DOUBLE( minRange ) ELSE minRange = FLOAT( minRange )
+ IF Keyword_Set(double) THEN maxRange = DOUBLE( maxRange ) ELSE maxRange = FLOAT( maxRange )
+
+ ; Make sure we have a valid range.
+ IF maxRange EQ minRange THEN Message, 'Range max and min are coincidental'
+
+ ; Check keyword parameters.
+ IF Keyword_Set(double) THEN BEGIN
+ IF N_Elements(vectorMin) EQ 0 THEN vectorMin = Double( Min(FPUFIX(vector), NAN=1) ) $
+ ELSE vectorMin = Double(vectorMin)
+ IF N_Elements(vectorMax) EQ 0 THEN vectorMax = DOUBLE( Max(FPUFIX(vector), NAN=1) ) $
+ ELSE vectorMax = DOUBLE( vectorMax )
+ ENDIF ELSE BEGIN
+ IF N_Elements(vectorMin) EQ 0 THEN vectorMin = FLOAT( Min(FPUFIX(vector), NAN=1) ) $
+ ELSE vectorMin = FLOAT( vectorMin )
+ IF N_Elements(vectorMax) EQ 0 THEN vectorMax = FLOAT( Max(FPUFIX(vector), NAN=Keyword_Set(nan)) ) $
+ ELSE vectorMax = FLOAT( vectorMax )
+ ENDELSE
+
+ ; Trim vector before scaling.
+ index = Where(Finite(vector) EQ 1, count)
+ IF count NE 0 THEN BEGIN
+ IF Keyword_Set(double) THEN trimVector = Double(vector) ELSE trimVector = Float(vector)
+ trimVector[index] = vectorMin > vector[index] < vectorMax
+ ENDIF ELSE BEGIN
+ IF Keyword_Set(double) THEN trimVector = vectorMin > Double(vector) < vectorMax ELSE $
+ trimVector = vectorMin > Float(vector) < vectorMax
+ ENDELSE
+
+ ; Calculate the scaling factors.
+ scaleFactor = [((minRange * vectorMax)-(maxRange * vectorMin)) / $
+ (vectorMax - vectorMin), (maxRange - minRange) / (vectorMax - vectorMin)]
+
+ ; Clear math errors.
+ void = Check_Math()
+
+ ; Return the scaled vector.
+ IF Keyword_Set(preserve_type) THEN BEGIN
+ RETURN, FPUFIX(Convert_To_Type(trimVector * scaleFactor[1] + scaleFactor[0], Size(vector, /TNAME)))
+ ENDIF ELSE BEGIN
+ RETURN, FPUFIX(trimVector * scaleFactor[1] + scaleFactor[0])
+ ENDELSE
+
+END
diff --git a/Patch_lib/patch_check_overlap.pro b/Patch_lib/patch_check_overlap.pro
new file mode 100644
index 0000000000000000000000000000000000000000..520d4d351f428bd09d8185d4c3768537ed745c76
--- /dev/null
+++ b/Patch_lib/patch_check_overlap.pro
@@ -0,0 +1,52 @@
+;last update: 2014-10-17 by Andrea La Camera
+;+
+;At this point, the overlap value depends on three different "sources":
+; (a) the one coming from the Enclosed Energy of the PSF
+; (b) the one coming from the difference (np - n) (i.e. size of the psf - size of the domain)
+; (c) the user-defined value
+;When called, the routine computes (a) and (b), checks (c) and
+;then verifies which value is the largest one.
+;Finally the routine prints warning messages if requested by the
+;keyword ask_confirm=1
+;-
+
+pro patch_check_overlap, par, parent_ID, auto_overlap, check_overlap, VERBOSE=verbose
+
+if KEYWORD_SET(verbose) then verb=1 else verb=0
+check_overlap = 0
+;(a) compute EE, find max(EE(psf_k))
+;
+;(b) find np and n
+ delta = 0
+ max_np = 0
+ if par.psf_type EQ 0 then psf=readfits((*par.psf_filename)[0],/SILENT)
+ for k=0,par.cols*par.rows-1 do begin
+ if par.psf_type EQ 1 then begin
+ psf_curr=readfits((*par.psf_filename)[k],/SILENT)
+ endif else psf_curr=psf[*,*,k]
+ NP=(size(psf_curr))[1]
+ MP=(size(psf_curr))[2]
+ delta_k=patch_psf_ee(psf_curr, par.threshold)
+ if (delta_k GT delta) then delta = delta_k
+ if (max([NP,MP]) GT max_np) then max_np=max([NP,MP])
+ endfor
+ ;max size of the domains
+ max_n = max([par.domx,par.domy])
+
+;(c) check par.overlap (the user-defined value)
+; and find the largest value
+ auto_overlap = max([delta, ROUND(((max_np-max_n)>0)/2)])
+ if (par.overlap LT auto_overlap) then check_overlap = 1
+
+;print warnings if necessary....
+ if (verb AND check_overlap) then begin
+ a=dialog_message(["THE OVERLAP IS TOO SMALL!", $
+ "The user-defined overlap must be greater than",$
+ "the auto value ("+strtrim(auto_overlap,1)+")",$
+ " ",$
+ "Select OK to continue and accept this number",$
+ "or Cancel to return to the GUI."],$
+ dialog_parent=parent_ID, /CANCEL, /INFO)
+ if (a EQ 'Cancel') then check_overlap=0
+ endif
+end
diff --git a/Patch_lib/patch_deconv.pro b/Patch_lib/patch_deconv.pro
new file mode 100644
index 0000000000000000000000000000000000000000..9462dc73b987deb6b841739fc62b22c6f63e3f26
--- /dev/null
+++ b/Patch_lib/patch_deconv.pro
@@ -0,0 +1,226 @@
+;last update: 2014-04-30 by Andrea La Camera
+;+
+;-
+
+pro patch_deconv, par, SAVE_OUT=save_out ;, obj
+
+ COMMON image, img, header, N, M, a_dim, Nori, Mori
+ COMMON background, bg ;bg must be N x M, as the input image "img" is.
+ COMMON output, out, res_array, residual ; reconstructed object with original dimension Nori x Mori
+ ; res_array is the residual statistics (one for each sub-domain)
+ ; residual is the final map of residuals
+ on_error, 2
+;save several information in header
+ fxaddpar, header, 'IMG ',par.img_filename,"Image: blurred image filename"
+ if par.bg_choice EQ 0 then begin
+ fxaddpar, header, 'SKY BG ',par.bg, "Image: sky background [ADU/pixel]"
+ endif else begin
+ fxaddpar, header, 'SKY BG ',par.bg_file, "Image: sky background from file"
+ endelse
+ fxaddpar, header, 'RON ',par.sigma_ron, "Image: RON value [e-/pixel]"
+ fxaddpar, header, 'GAIN ',par.gain, "Image: GAIN value [e-/ADU]"
+ fxaddpar, header, 'N_COLS ',par.cols, "Sub-domain: number of columns"
+ fxaddpar, header, 'N_ROWS ',par.rows, "Sub-domain: number of rows"
+ fxaddpar, header, 'OVERLAP ',par.overlap, "Sub-domain: overlap"
+ fxaddpar, header, 'PSF_TYPE',par.psf_type, "PSF: type (0=cube, 1=single files)"
+ if par.psf_type EQ 0 then begin
+ fxaddpar, header, 'PSF_CUBE',(*par.psf_filename)[0], "PSF: filename"
+ endif else begin
+ for k=0,par.cols*par.rows-1 do begin
+ fxaddpar, header, 'PSF#'+strtrim(k,1),(*par.psf_filename)[k], "PSF: filename"
+ endfor
+ endelse
+ fxaddpar, header, 'EE_THRES',par.threshold, "PSF: Enclosed Energy threshold"
+ fxaddpar, header, 'DECONV ',par.method, "Deconvolution: method (0=RL,1=SGP)"
+ fxaddpar, header, 'STOP_CRI',par.stop_crit, "Deconvolution: stopping rule"
+ fxaddpar, header, 'TOL ',par.tol, "Deconvolution: Tolerance (if requested)"
+
+;evaluate the median of the background
+ bg_med=median(bg)
+
+;zero-padding of image (in order to consider overlap)
+ NN=N+2*par.overlap
+ MM=M+2*par.overlap
+ img2=zero_padding(img,NN,MM)
+;obj2=zero_padding(obj,NN,MM)
+;zero-padding of the background array
+ bg2=zero_padding(bg,NN,MM)
+
+;add the median bg value to the zero-padded img
+ index=where(img2 LE 0, count)
+ if count GT 0 then img2[index]=bg_med
+
+;add the median bg value to the zero-padded bg
+ index=where(bg2 LE 0, count)
+ if count GT 0 then bg2[index]=bg_med
+
+;add the Read-Out Noise variance to img2 and bg2 (Gaussian+Poisson noise correction)
+ img2+=((par.sigma_ron/par.gain)^2.)
+ bg2+=((par.sigma_ron/par.gain)^2.)
+ bg_med+=((par.sigma_ron/par.gain)^2.)
+
+
+;output image (the reconstructed object), residual and other stuff
+ out=make_array(N,M,/float)
+ residual=make_array(N,M,/float)
+ res_array=make_array(par.cols*par.rows,5)
+
+;dimension of each sub-domain (+ overlap)
+ NS=par.domx+2*par.overlap
+ MS=par.domy+2*par.overlap
+
+;read psfs
+ if par.psf_type EQ 0 then psf=readfits((*par.psf_filename)[0],/SILENT)
+
+;read MASK for PROJECTION (in SGP)
+ if (par.mask_choice) then begin
+ pr_mask=readfits(par.mask_filename, /SILENT)
+ pr_mask=zero_padding(pr_mask,NN,MM,value=1.)
+ endif else begin
+ pr_mask=make_array(NN,MM,value=1.)
+ endelse
+
+ tt=systime(1)
+;-------> RUN!
+ for k=0,par.cols*par.rows-1 do begin
+ t0=systime(1)
+ ;cycle on sub-domains
+ ;selection of current sub-domain
+ x0=(k mod par.cols)*par.domx
+ y0=(k / par.cols)*par.domy
+ x1=x0+NS-1
+ y1=y0+MS-1
+ g=img2[x0:x1,y0:y1]
+ b=bg2[x0:x1,y0:y1]
+ sub_pr_mask=pr_mask[x0:x1,y0:y1]
+ ;o=obj2[x0:x1,y0:y1]
+ ;flux evaluation
+ flux=total(g-b, /DOUBLE)
+ ;selection of current PSF
+ if par.psf_type EQ 1 then begin
+ psf_curr=readfits((*par.psf_filename)[k],/SILENT)
+ endif else psf_curr=psf[*,*,k]
+ ;(re-)normalization of the psf
+ psf_curr=psf_curr/total(psf_curr,/double)
+ NP=(size(psf_curr))[1]
+ MP=(size(psf_curr))[2]
+ ;compute Delta (the Enclosed Energy crosses the threshold at distance = Delta)
+;delta=patch_psf_ee(psf_curr, par.threshold)
+ ;print, "DELTA = ", delta
+;NS1=NS+2*delta
+;MS1=MS+2*delta
+ case par.method of
+ 0: begin
+ ;RICHARDSON-LUCY METHOD + BOUNDARY
+ ;resizing
+;dim=max([NS1,MS1,NP,MP])
+ dim=max([NS,MS,NP,MP])
+ dim= dim + dim mod 2 ; I want an even number !
+ if dim LT 2l^ceil(alog(dim)/alog(2.)) then dim = 2l^ceil(alog(dim)/alog(2.))
+ psf_curr=zero_padding(psf_curr, dim, dim)
+ g=zero_padding(g,dim,dim)
+ b=zero_padding(b,dim,dim,value=bg_med)
+ sub_pr_mask=zero_padding(sub_pr_mask,dim,dim,value=1.)
+ print, "Working dimensions = ", dim
+ ;FFT of the PSF
+ tf_psf=FFT(shift(psf_curr,dim/2,dim/2),-1,/double)*dim*dim
+ ;computation of the initial arrays
+ mask=make_array(x1-x0+1,y1-y0+1, value=1.)
+ mask=zero_padding(mask,dim,dim)
+ alpha=FLOAT(FFT(conj(tf_psf)*FFT(mask,-1,/double),+1,/double))
+ w=1./alpha
+ index=where(alpha LT 1e-2, count)
+ if count GT 0 then w[index]=0.
+ ;initial rec --- with exact flux in the reconstruction region
+ rec=make_array(dim,dim,value=(flux/total(alpha,/DOU)))
+ if count GT 0 then rec[index]=0.
+ rec=rec*sub_pr_mask
+ ;run RL iterations
+ iter=0
+ print, 'DECONV RL --> nb of iterations = ', (*par.tab_iter)[k]
+ patch_deconv_rl, g,tf_psf,b,w,(*par.tab_iter)[k],par.stop_crit,rec, iter, res, res_stat
+ ;save rec'd object into correct
+ ;position
+ out[x0:x0+par.domx-1,y0:y0+par.domy-1]=rec[(dim-par.domx)/2:(dim+par.domx)/2-1,$
+ (dim-par.domy)/2:(dim+par.domy)/2-1]
+ ; save residual into array
+ residual[x0:x0+par.domx-1,y0:y0+par.domy-1]=res[(dim-par.domx)/2:(dim+par.domx)/2-1,$
+ (dim-par.domy)/2:(dim+par.domy)/2-1]
+ res_array[k,*]=res_stat
+ print, " --> sub-domain #"+strtrim(k,1)+" done! (elapsed time:"+$
+ strtrim(systime(1)-t0,1)+" sec) (it = "+strtrim(iter,1)+")"
+ end
+ 1: begin
+ ; SGP METHOD
+ ; if g is very small (i.e. less than 2^(N-1)) while the PSF is less than
+ ; or equal to 2^N) . For istance: g=200x200 and PSF = 400x400 -->
+ ; g=256x256 and PSF = 512x512. Next, SGP will extend g to the actual size of
+ ; the PSF (512).
+;if (max([NS1,MS1]) LE 2^floor(alog(max([NP,MP]))/alog(2))) then begin
+ if (max([NS,MS]) LE 2^floor(alog(max([NP,MP]))/alog(2))) then begin
+ g=zero_padding(g,2^floor(alog(max([NP,MP]))/alog(2)),2^floor(alog(max([NP,MP]))/alog(2)),value=bg_med)
+ b=zero_padding(b,2^floor(alog(max([NP,MP]))/alog(2)),2^floor(alog(max([NP,MP]))/alog(2)),value=bg_med)
+ sub_pr_mask=zero_padding(sub_pr_mask,2^floor(alog(max([NP,MP]))/alog(2)),2^floor(alog(max([NP,MP]))/alog(2)),value=1.)
+ ;o=zero_padding(o,2^floor(alog(max([NP,MP]))/alog(2)), 2^floor(alog(max([NP,MP]))/alog(2)))
+; NS1=(size(g))[1]
+; MS1=(size(g))[2]
+ NS=(size(g))[1]
+ MS=(size(g))[2]
+; dim=max([NS1,MS1,NP,MP])
+ dim=max([NS,MS,NP,MP])
+ if dim LT 2l^ceil(alog(dim)/alog(2.)) then dim = 2l^ceil(alog(dim)/alog(2.))
+ psf_curr=zero_padding(psf_curr, dim, dim)
+ ;o=zero_padding(o, dim, dim)
+ endif else begin ; in all other cases enlarge g,b,and PSF to the same dim.
+; dim=max([NS1,MS1,NP,MP])
+ dim=max([NS,MS,NP,MP])
+ dim= dim + dim mod 2 ; I want an even number !
+ psf_curr=zero_padding(psf_curr, dim, dim)
+ g=zero_padding(g,dim,dim,value=bg_med)
+ b=zero_padding(b,dim,dim,value=bg_med)
+ sub_pr_mask=zero_padding(sub_pr_mask,dim,dim,value=1.)
+ ;o=zero_padding(o,dim,dim)
+ endelse
+ print, "Working dimensions = ", dim
+ iter=0
+ ;print, 'DECONV SGP --> nb of iterations = ', (*par.tab_iter)[k]
+ patch_deconv_sgp, g, psf_curr, b, (*par.tab_iter)[k], $
+ par.stop_crit, par.tol, rec, iter, sub_pr_mask, res, res_stat
+ ;save rec'd object into correct
+ ;position
+ iter=iter+1
+ dim=(size(rec))[1]
+ out[x0:x0+par.domx-1,y0:y0+par.domy-1]=rec[(dim-par.domx)/2:(dim+par.domx)/2-1,$
+ (dim-par.domy)/2:(dim+par.domy)/2-1]
+ ; save residual into array
+ residual[x0:x0+par.domx-1,y0:y0+par.domy-1]=res[(dim-par.domx)/2:(dim+par.domx)/2-1,$
+ (dim-par.domy)/2:(dim+par.domy)/2-1]
+
+ res_array[k,*]=res_stat
+ print, " --> sub-domain #"+strtrim(k,1)+" done! (elapsed time:"+$
+ strtrim(systime(1)-t0,1)+" sec) (it = "+strtrim(iter-1,1)+")"
+ ;print, "err = ", err
+ ;print, "Discr = ", Discr
+ end
+ else: begin
+ message, "Deconvolution method not implemented!"
+ endelse
+ endcase
+;add other info to header
+ fxaddpar, header, "DOM#"+strtrim(k,1),iter-1 , "Deconvolution: iterations in "+strtrim(systime(1)-t0,1)+"(sec)"
+ endfor
+;last header update
+ fxaddpar,header, "DEC_TIME",strtrim(systime(1)-tt,1),"Deconvolution: Total time (sec)"
+
+;crop reconstructed image to input (original) dimensions
+ out=out[(N-Nori)/2:(N+Nori)/2-1, (M-Mori)/2:(M+Mori)/2-1]
+;crop the residual image to input (original) dimensions
+ residual=residual[(N-Nori)/2:(N+Nori)/2-1, (M-Mori)/2:(M+Mori)/2-1]
+
+;save output in FITS file if SAVE_OUT keyword is set
+ if KEYWORD_SET(save_out) then begin
+ writefits, par.out_filename, out, header
+ writefits, par.res_filename, residual, header
+ endif
+;ok! done!
+end
diff --git a/Patch_lib/patch_deconv_rl.pro b/Patch_lib/patch_deconv_rl.pro
new file mode 100644
index 0000000000000000000000000000000000000000..6a87e28f443c0ac2b07e29864bfde87ce06c38f9
--- /dev/null
+++ b/Patch_lib/patch_deconv_rl.pro
@@ -0,0 +1,16 @@
+pro patch_deconv_rl, g,tf_psf, bg, w, tot_iter,stop_crit, rec, iter, res, res_stat
+
+while (iter LT tot_iter) AND (stop_crit EQ 0) do begin
+ dummy=g/(float(FFT(FFT(rec,-1,/double)*tf_psf,+1,/double))+bg)
+ dummy=float(FFT(conj(tf_psf)*FFT(dummy,-1,/double),+1,/double))
+ rec=rec*dummy*w
+ iter++
+ ;print, "Iteration # ", iter
+endwhile
+;compute the residual array:
+;R_j = [g_j - (A_j*f_j + b_j)] / sqrt (A_j*f_j + b_j)
+dummy=(float(FFT(FFT(rec,-1,/double)*tf_psf,+1,/double))+bg)
+res=(g-dummy)/sqrt(dummy)
+histogauss, res, res_stat, /NOPLOT
+
+end
diff --git a/Patch_lib/patch_deconv_sgp.pro b/Patch_lib/patch_deconv_sgp.pro
new file mode 100644
index 0000000000000000000000000000000000000000..67199d6e0f9a0a692956cb8eee7a3567ab01f4a8
--- /dev/null
+++ b/Patch_lib/patch_deconv_sgp.pro
@@ -0,0 +1,13 @@
+pro patch_deconv_sgp, gn, A, bg, maxit, stopcrit, tol, x, iter, pr_mask, res, res_stat
+
+;call sgp_deblurring.pro and Afunction.pro (part of PRISMA Software)
+;see the header of the files for more details and license.
+
+boundary=1 ;YES!
+verbose = 0
+sgp_deblurring, A, gn, x, iter, err, Discr, times, pr_mask, res, res_stat, $
+ bg = bg, bound=boundary, $
+ maxit = maxit, stopcrit = stopcrit+1, tol = tol, $
+ verb = verbose
+
+end
diff --git a/Patch_lib/patch_gui.pro b/Patch_lib/patch_gui.pro
new file mode 100644
index 0000000000000000000000000000000000000000..49b1033a52ccdb493ccf45057fccaf10b3af5292
--- /dev/null
+++ b/Patch_lib/patch_gui.pro
@@ -0,0 +1,710 @@
+;PATCH_GUI.PRO
+;This routine manages the Graphical User Interface (GUI)
+;and its events for the software PATCH
+;Andrea La Camera (DIBRIS, Univ. of Genova), November 2011
+;last modified: 2014-10-17
+;---------
+;
+
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+pro patch_gui_event, event
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+ COMMON image, img, header, N, M, a_dim, Nori, Mori
+ COMMON background, bg
+ COMMON working_dir, wd
+ COMMON output, out, res_array, residual ; reconstructed object with original dimension Nori x Mori
+
+
+ WIDGET_CONTROL, event.ID, GET_UVALUE=choice
+ WIDGET_CONTROL, event.TOP, GET_UVALUE=state
+; handle a kill request (considered as a cancel event).
+ if tag_names(event, /STRUCTURE_NAME) eq 'WIDGET_KILL_REQUEST' then begin
+ widget_control, event.top, /DESTROY
+ endif
+
+ case choice of
+ 'img_load' : begin
+ state.par.img_filename=dialog_pickfile(dialog_parent=event.top, $
+ FILTER=['*.fits', '*.fit', '*.fts', '*.*'], $
+ PATH=wd, GET_PATH=new_wd, $
+ /MUST_EXIST, title="Load a FITS file (must be a 2d array)")
+ if state.par.img_filename EQ '' then return
+ wd=new_wd
+ a=patch_load_image(state, event.top, /RESET, /VERBOSE)
+ if a then patch_show_image, state, /redraw else return
+ end
+ 'img_define' : begin
+ if state.par.img_filename EQ '' then return
+ result=patch_img_define("SUB-DOMAINS DEFINITION", state, group_leader=event.top)
+ state.par=result
+ if (N GT (size(img))[1]) OR (M GT (size(img))[2]) then begin
+ img = zero_padding(img, N, M)
+ endif
+ info=["filename: "+state.par.img_filename, $
+ "size: "+strtrim(N,1)+' x '+strtrim(M,1), $
+ "domains ("+strtrim(state.par.rows*state.par.cols,1)+"): size "+$
+ strtrim(state.par.domx,1)+' x '+strtrim(state.par.domy,1), $
+ "current domain: none"]
+ ;show info
+ WIDGET_CONTROL, state.id.img_info, SET_VALUE=info
+ patch_show_image, state, /REDRAW
+ end
+ 'iarea' : begin
+ if state.par.img_filename EQ '' then return
+ ;print, event.x, event.y, event.clicks, event.press, event.modifiers
+ if (event.press EQ 1) AND (event.modifiers EQ 0) then begin
+ ;select a sub-domain
+ if state.par.rows*state.par.cols GT 1 then begin
+ patch_show_image, state ;Clean area
+ x0=(FIX(event.x/state.par.img_zoom-state.par.overlap)>0<(N-1)) / state.par.domx
+ y0=(FIX(event.y/state.par.img_zoom-state.par.overlap)>0<(M-1)) / state.par.domy
+ ;draw the box around selected sub-region
+ plots, (indgen(state.par.domx+2*state.par.overlap)+x0*state.par.domx)*state.par.img_zoom, $
+ (y0*state.par.domy)*state.par.img_zoom, /DEV, $
+ color='FFFFFF'x, LINESTYLE=0, thick=2
+ plots, (indgen(state.par.domx+2*state.par.overlap)+x0*state.par.domx)*state.par.img_zoom, $
+ ((y0+1)*state.par.domy+2*state.par.overlap-2)*state.par.img_zoom, /DEV, $
+ color='FFFFFF'x, LINESTYLE=0, thick=2
+ plots, (x0*state.par.domx)*state.par.img_zoom, $
+ (indgen(state.par.domy+2*state.par.overlap)+y0*state.par.domy)*state.par.img_zoom, $
+ /DEV, color='FFFFFF'x, LINESTYLE=0, thick=2
+ plots, ((x0+1)*state.par.domx+2*state.par.overlap-2)*state.par.img_zoom, $
+ (indgen(state.par.domy+2*state.par.overlap)+y0*state.par.domy)*state.par.img_zoom, $
+ /DEV, color='FFFFFF'x, LINESTYLE=0, thick=2
+ if ((state.par.domx+state.par.domy)*state.par.img_zoom) LT 51 then begin
+ xyouts, ((x0+0.5)*state.par.domx+state.par.overlap)*state.par.img_zoom, $
+ ((y0+0.47)*state.par.domy+state.par.overlap)*state.par.img_zoom, $
+ strtrim(y0*state.par.cols+x0,1), $
+ /DEVICE, charsize=1.2, ALIGNMENT=0.5, charthick=2.
+ endif
+ ;show info
+ info=["filename: "+state.par.img_filename, $
+ "size: "+strtrim(N,1)+' x '+strtrim(M,1), $
+ "domains ("+strtrim(state.par.rows*state.par.cols,1)+"): size "+$
+ strtrim(state.par.domx,1)+' x '+strtrim(state.par.domy,1), $
+ "current domain: "+strtrim(y0*state.par.cols+x0,1)]
+ WIDGET_CONTROL, state.id.img_info, SET_VALUE=info
+
+ endif
+ endif
+ if (event.press EQ 4) then begin
+ ;select options (zoom, scale, color, etc)
+ result=patch_img_options("OPTIONS", state, group_leader=event.top)
+ if result.img_zoom EQ state.par.img_zoom then begin
+ state.par.img_scale = result.img_scale
+ state.par.img_color = result.img_color
+ patch_show_image, state ; REDRAW is not necessary
+ endif else begin
+ state.par.img_zoom = result.img_zoom
+ state.par.img_scale = result.img_scale
+ state.par.img_color = result.img_color
+ patch_show_image, state, /REDRAW ; REDRAW is necessary
+ endelse
+ endif
+ end
+ 'psf_type': begin
+ state.par.psf_type=event.value
+ ;reset info of psf
+ widget_control, state.id.psf_info, set_value=strarr(4)
+ patch_set_gui, state
+ end
+ 'psf_load' : begin
+ nb_of_psfs = state.par.cols*state.par.rows
+ if state.par.psf_type EQ 0 then begin
+ res=dialog_pickfile(FILTER=['*.fits', '*.fit', '*.fts', '*.*'], /MUST_EXIST, $
+ path=wd, get_path=new_wd, $
+ title="Load the FITS file containing the "+strtrim(nb_of_psfs,1)+" PSFs.", $
+ dialog_parent=event.top)
+ if res EQ '' then return
+ wd=new_wd
+ (*state.par.psf_filename)[0]=res
+ a=patch_load_psfs(state, event.top)
+ endif
+ if state.par.psf_type EQ 1 then begin
+ res=dialog_pickfile(FILTER=['*.fits', '*.fit', '*.fts', '*.*'], $
+ path=wd, get_path=new_wd, $
+ /MUST_EXIST, /MULTIPLE_FILES, $
+ title="Load "+strtrim(nb_of_psfs,1)+$
+ " FITS files (must be 2d arrays)", dialog_parent=event.top)
+ if (n_elements(res) EQ 0) OR (res[0] EQ '') then return
+ if n_elements(res) NE nb_of_psfs then begin
+ a=dialog_message("Wrong number of PSFs: "+strtrim(nb_of_psfs,1)+$
+ " files expected.", dialog_parent=event.top, /ERROR)
+ return
+ endif
+ wd=new_wd
+ ;load psfs
+ (*state.par.psf_filename)=res
+ a=patch_load_psfs(state, event.top)
+ endif
+ end
+ 'threshold': state.par.threshold = float(event.value)/100.
+ 'auto_overlap': begin
+ if state.par.img_filename EQ '' then return
+ print, "AUTO OVERLAP calculating ..."
+ patch_check_overlap, state.par, event.top, auto_overlap, check_overlap, VERBOSE=0
+ state.par.overlap=auto_overlap
+ WIDGET_CONTROL, state.id.overlap, SET_VALUE=state.par.overlap
+ print, "...done!"
+ end
+ 'overlap' : begin
+ state.par.overlap=event.value
+ print, state.par.overlap
+ end
+ 'check_overlap' : begin
+ if state.par.img_filename EQ '' then return
+ print, "CHECK OVERLAP calculating ... "
+ patch_check_overlap, state.par, event.top, auto_overlap, check_overlap, /VERBOSE
+ if check_overlap then begin
+ state.par.overlap=auto_overlap
+ WIDGET_CONTROL, state.id.overlap, SET_VALUE=state.par.overlap
+ endif
+ print, "...done!"
+ end
+ 'redraw' : patch_show_image, state, /REDRAW
+ 'bg_choice': begin
+ state.par.bg_choice = event.value
+ patch_set_gui, state
+ end
+ 'bg' : begin
+ state.par.bg=event.value
+ bg = make_array(N,M,value=state.par.bg)
+ end
+ 'bg_file' : begin
+ state.par.bg_file=dialog_pickfile(dialog_parent=event.top, $
+ FILTER=['*.fits', '*.fit', '*.fts', '*.*'], $
+ PATH=wd, GET_PATH=new_wd, $
+ /MUST_EXIST, title="Load a FITS file (must be a 2d array)")
+ if state.par.bg_file EQ '' then return
+ a=patch_load_background(state, event.top, /VERB)
+ wd=new_wd
+ end
+ 'sigma_ron': state.par.sigma_ron = event.value
+ 'gain' : state.par.gain = event.value
+ 'method' : begin
+ state.par.method=event.index
+ patch_set_gui, state
+ end
+ 'stop_crit' : begin
+ state.par.stop_crit = event.index
+ patch_set_gui, state
+ end
+ 'tol' : state.par.tol = event.value
+ 'mask_choice': begin
+ state.par.mask_choice = event.select
+ patch_set_gui, state
+ end
+ 'mask_filename': begin
+ state.par.mask_filename=dialog_pickfile(dialog_parent=event.top, $
+ FILTER=['*.fits', '*.fit', '*.fts', '*.*'], $
+ PATH=wd, GET_PATH=new_wd, $
+ /MUST_EXIST, title="Load a FITS file (must be a 2d array)")
+ if state.par.mask_filename EQ '' then return
+ wd=new_wd
+ mask=readfits(state.par.mask_filename, /SILENT)
+ dim_mask_1=(size(mask))[1]
+ info=["filename: "+state.par.mask_filename, $
+ "size: "+strtrim((size(mask))[1],1)+' x '+strtrim((size(mask))[2],1)]
+ ;show info
+ WIDGET_CONTROL, state.id.mask_info, SET_VALUE=info
+ end
+ 'it_choice': begin
+ state.par.it_choice=event.value
+ patch_set_gui, state
+ ;compute the NoI table here (1)
+ patch_noi_table_compute, event.value, state.par
+ end
+ 'tot_iter': begin
+ state.par.tot_iter=event.value
+ ;compute the NoI table also here (2)
+ patch_noi_table_compute, 0, state.par
+ end
+ 'it4': begin
+ widget_control, state.id.it4, get_value=dummy
+ if dummy[event.x,event.y] LE 0 then begin
+ a=dialog_message("Insert positive (>0) values!", /ERROR)
+ endif
+ dummy=float(rotate(dummy,7))
+ state.par.it4=dummy
+ ;compute the NoI table also here (3)
+ patch_noi_table_compute, 1, state.par
+ end
+ 'itedit':begin
+ if n_elements(*state.par.tab_iter) NE state.par.cols*state.par.rows then begin
+ ;compute the NoI table also here (4)
+ patch_noi_table_compute, state.par.it_choice, state.par
+ endif
+ res=patch_noi_table(state.par,GROUP_LEADER=event.top)
+ (*state.par.tab_iter)=(*res.tab_iter)
+ ;print, '--> EXIT : ', *state.par.tab_iter
+ end
+ 'go' : begin
+ if n_elements(img) EQ 0 then begin
+ a=dialog_message("Input image undefined!", /ERROR)
+ return
+ endif
+ if n_elements(bg) EQ 0 then begin
+ a=dialog_message("Background undefined! Insert positive (>0) values", /ERROR)
+ return
+ endif
+ if (*state.par.psf_filename)[0] EQ '' then return
+ ;;;;; MODIFICA PER RMSE in ogni sotto-dominio
+ ;a= dialog_pickfile(dialog_parent=event.top, FILTER=['*.fits'], $
+ ; path=wd, get_path=new_wd, $
+ ; /MUST_EXIST, title="CARICA L'OGGETTO NORMALIZZATO")
+ ;print, a
+ ;if a NE '' then begin
+ ; fits_read, a, obj
+ ; N_tmp=(size(obj))[1]
+ ; M_tmp=(size(obj))[2]
+ ; ;check if zero_padding is neccessary...
+ ; if (state.par.cols*state.par.domx GT N_tmp) OR $
+ ; (state.par.rows*state.par.domy GT M_tmp) then begin
+ ; obj = zero_padding(obj, N, M)
+ ; endif
+ ;endif else obj=make_array(N,M)
+ ;help, obj
+ ;;;;;;;fine modifica
+ t0=systime(1)
+ print, "Please, wait..."
+ WIDGET_CONTROL, state.id.run_deconv, SET_VALUE=" .... wait .... "
+ WIDGET_CONTROL, state.id.run_deconv, SENSITIVE=0
+ patch_deconv, state.par ;, obj ;;;; aggiunto obj !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
+ WIDGET_CONTROL, state.id.run_deconv, SET_VALUE=" RUN DECONVOLUTION "
+ WIDGET_CONTROL, state.id.run_deconv, /SENSITIVE
+ print, "...done!"
+ print, "Total deconvolution time (s) = ", systime(1)-t0
+ print, "Average time for each sub-domain = ", $
+ (systime(1)-t0)/(state.par.rows*state.par.cols)
+ patch_show_out, state, /REDRAW ; REDRAW is necessary
+ end
+ 'save_out' : begin
+ res=dialog_pickfile(dialog_parent=event.top, $
+ FILTER=['*.fits', '*.fit', '*.fts'], $
+ path=wd, get_path=new_wd, $
+ /OVERWRITE_PROMPT, /WRITE, $
+ title="Save as a FITS file (must be a 2d array)")
+ if res EQ '' then return else state.par.out_filename=res
+ wd=new_wd
+ if n_elements(Nori) EQ 0 then begin
+ info=["filename: "+state.par.out_filename, $
+ "size: not yet defined "]
+ endif else begin
+ info=["filename: "+state.par.out_filename, $
+ " size: "+strtrim(Nori,1)+' x '+strtrim(Mori,1)]
+ endelse
+ ;show info (OK or ERROR)
+ WIDGET_CONTROL, state.id.out_info, SET_VALUE=info
+ if n_elements(out) GT 0 then begin
+ ;save output in FITS file
+ writefits, state.par.out_filename, out, header
+ endif
+ end
+ 'show_out': begin
+ if n_elements(out) GT 0 then patch_show_out, state, /REDRAW ; REDRAW is necessary
+ end
+ 'header': begin
+ if n_elements(header) GT 0 then xdispstr, header, TITLE="FITS HEADER"
+ end
+ 'save_res' : begin
+ res=dialog_pickfile(dialog_parent=event.top, $
+ FILTER=['*.fits', '*.fit', '*.fts'], $
+ path=wd, get_path=new_wd, $
+ /OVERWRITE_PROMPT, /WRITE, $
+ title="Save as a FITS file (must be a 2d array)")
+ if res EQ '' then return else state.par.res_filename=res
+ wd=new_wd
+ if n_elements(Nori) EQ 0 then begin
+ info=["filename: "+state.par.res_filename, $
+ "size: not yet defined "]
+ endif else begin
+ info=["filename: "+state.par.res_filename, $
+ " size: "+strtrim(Nori,1)+' x '+strtrim(Mori,1)]
+ endelse
+ ;show info (OK or ERROR)
+ WIDGET_CONTROL, state.id.res_info, SET_VALUE=info
+ if n_elements(residual) GT 0 then begin
+ ;save output in FITS file
+ writefits, state.par.res_filename, residual, header
+ endif
+ end
+ 'show_res': begin
+ if n_elements(residual) GT 0 then patch_show_res, state, /REDRAW ; REDRAW is necessary
+ end
+ 'stat_res': begin
+ if n_elements(residual) GT 0 then begin
+ window,/free,title="Residual Statistics - HISTOGAUSS RESULT"
+ histogauss, residual, stat, charsize=1.8
+ endif
+ end
+ 'oarea': begin
+ if state.par.img_filename EQ '' then return
+ if (event.press EQ 4) then begin
+ ;select options (zoom, scale, color, etc)
+ result=patch_img_options("OPTIONS", state, group_leader=event.top)
+ if result.img_zoom EQ state.par.img_zoom then begin
+ state.par.img_scale = result.img_scale
+ state.par.img_color = result.img_color
+ patch_show_out, state ; REDRAW is not necessary
+ endif else begin
+ state.par.img_zoom = result.img_zoom
+ state.par.img_scale = result.img_scale
+ state.par.img_color = result.img_color
+ patch_show_out, state, /REDRAW ; REDRAW is necessary
+ endelse
+ endif
+ end
+ 'load_par': begin
+ a=dialog_pickfile(dialog_parent=event.top, FILTER=['*.sav'], $
+ path=wd, get_path=new_wd, $
+ /MUST_EXIST, /READ)
+ if a EQ '' then return
+ wd=new_wd
+ restore, a, /VERBOSE
+ ;restore all parameters and set the gui::: Read structure is "PARAM"
+ state.par=param
+ ; load image and define img, N, M
+ a=patch_load_image(state, event.top, /VERBOSE)
+ if a EQ 0 then return
+ ;check if zero_padding is neccessary...
+ if (state.par.cols*state.par.domx GT N) OR $
+ (state.par.rows*state.par.domy GT M) then begin
+ N=state.par.cols*state.par.domx
+ M=state.par.rows*state.par.domy
+ img = zero_padding(img, N, M)
+ endif
+ ;show image
+ patch_show_image, state, /REDRAW
+ ;info about image
+ info=["filename: "+state.par.img_filename, $
+ "size: "+strtrim(N,1)+' x '+strtrim(M,1), $
+ "domains ("+strtrim(state.par.rows*state.par.cols,1)+"): size "+$
+ strtrim(state.par.domx,1)+' x '+strtrim(state.par.domy,1), $
+ "current domain: none"]
+ ;show info
+ WIDGET_CONTROL, state.id.img_info, SET_VALUE=info
+ ; load PSF
+ WIDGET_CONTROL, state.id.psf_type, SET_VALUE=state.par.psf_type
+ a=patch_load_psfs(state, event.top)
+ WIDGET_CONTROL, state.id.threshold, SET_VALUE=fix(state.par.threshold*100.)
+ WIDGET_CONTROL, state.id.overlap, SET_VALUE=state.par.overlap
+ ;set "DECONVOLUTION" tab
+ if state.par.bg_choice EQ 0 then begin
+ WIDGET_CONTROL, state.id.bg, SET_VALUE=state.par.bg
+ bg = make_array(N,M,value=state.par.bg)
+ endif else begin
+ a=patch_load_background(state, event.top, /VERBOSE)
+ endelse
+ WIDGET_CONTROL, state.id.bg_choice, SET_VALUE=state.par.bg_choice
+ WIDGET_CONTROL, state.id.sigma_ron, SET_VALUE=state.par.sigma_ron
+ WIDGET_CONTROL, state.id.gain, SET_VALUE=state.par.gain
+ WIDGET_CONTROL, state.id.method, SET_COMBOBOX_SELECT=state.par.method
+ WIDGET_CONTROL, state.id.stop_crit, SET_COMBOBOX_SELECT=state.par.stop_crit
+ WIDGET_CONTROL, state.id.tol, SET_VALUE=state.par.tol
+
+
+
+ WIDGET_CONTROL, state.id.tot_iter, SET_VALUE=state.par.tot_iter
+ info=["filename: "+state.par.out_filename, $
+ " size: "+strtrim(Nori,1)+' x '+strtrim(Mori,1)]
+ WIDGET_CONTROL, state.id.out_info, SET_VALUE=info
+ info=["filename: "+state.par.res_filename, $
+ " size: "+strtrim(Nori,1)+' x '+strtrim(Mori,1)]
+ WIDGET_CONTROL, state.id.res_info, SET_VALUE=info
+ patch_set_gui, state
+ end
+; 'load_rec': begin
+; print, event.value
+; end
+ 'save_par': begin
+ a=dialog_pickfile(dialog_parent=event.top, FILTER=['*.sav'], $
+ path=wd, get_path=new_wd, $
+ /WRITE)
+ if a EQ '' then return
+ wd=new_wd
+ param=state.par
+ save, param, FILENAME=a, /VERBOSE
+ end
+; 'help' : begin
+; a=dialog_message("NOT YET IMPLEMENTED", /INFO, dialog_parent=event.top)
+ ;XDISPLAYFILE, nomefile, DONE_BUTTON="Done", Title="PATCH Help"
+; end
+ 'web' : begin
+ online_help, book="website.html"
+ end
+ 'about' : a=dialog_message([" PATCH ", $
+ " version 0.7 ", $
+ " written by Andrea La Camera ", $
+ " (DIBRIS, University of Genova) ", $
+ " February 2014 "], $
+ /INFO, title="About this software", $
+ dialog_parent=event.top)
+ 'exit' : begin
+ widget_control, event.top, /DESTROY
+ return
+ end
+
+ else: return
+ endcase ;end of main CASE
+
+; write the GUI state structure
+ WIDGET_CONTROL, event.top, SET_UVALUE=state
+end
+
+
+;;;;;;;;;;;;;;;;;;
+pro patch_gui
+;;;;;;;;;;;;;;;;;;
+
+ COMMON image, img, header, N, M, a_dim
+ COMMON working_dir, wd
+ COMMON output, out, res_array, residual ; reconstructed object with original dimension Nori x Mori
+
+ cd, current=wd
+ device, decompose=0
+
+;Definition of the structure "STATE", made of two sub-structures: "ID" and "PAR"
+
+ id = { img_define:0L, img_info:0L, $
+ ibase:0L, ikill:0L, iarea:0L, $
+ psf_type : 0L, psf_info:0L, text:0L, threshold:0L, $
+ overlap:0L, $
+ method:0L, stop_crit:0L, tol:0L, tot_iter:0L, $
+ mask_choice:0L, mask_filename:0L, mask_info:0L, $
+ bg:0L, bg_choice:0L, bg_file:0L, bg_info:0L, $
+ it_choice:0L, it4:0L, itedit:0L, $
+ sigma_ron:0L, gain:0L, run_deconv:0L, $
+ obase:0L, okill:0L, oarea:0L, $
+ out_info:0L, res_info:0L }
+
+ par = { img_filename:"", img_scale:0, img_zoom:1., $
+ img_color:0, psf_type:0, psf_filename:PTR_NEW(strarr(4)), $
+ threshold:0.95, $
+ domx:-1, domy:-1, rows:1, cols:1, overlap:0, $
+ method:0, stop_crit:0, tol:1e-4, tot_iter:1L, out_filename:"", $
+ mask_choice:0, mask_filename:"", $
+ bg:0., bg_file:'', bg_choice: 0B, $
+ it_choice:0, it4:[[1L,1L],[1L,1L]], tab_iter:PTR_NEW(1), $
+ sigma_ron:10., gain:1., res_filename:"" }
+
+ state={id:id, par:par}
+
+;Graphical User Interface starts here
+ base_id = widget_base(TITLE="PATCH - Graphical User Interface", MODAL=0, /COL, $
+ GROUP_LEADER=group, MBAR=bar)
+ menu1=WIDGET_BUTTON(bar, value="File", /MENU)
+ menu1b10=WIDGET_BUTTON(menu1, value="Load parameters...", uvalue="load_par")
+;OPENR, lun, file_which('_filelist.txt'), /GET_LUN
+;load_rec_desc = '1\Load recent '
+;line = ''
+;WHILE NOT EOF(lun) DO BEGIN
+; READF, lun, line
+; load_rec_desc = [load_rec_desc, '0\'+line]
+;ENDWHILE
+;close,lun
+;FREE_LUN, lun
+;load_rec_desc=[load_rec_desc,'2\Clear menu']
+;menu1b11=CW_PDMENU(menu1, load_rec_desc, /RETURN_NAME, uvalue="load_rec", /MBAR)
+ menu1b20=WIDGET_BUTTON(menu1, value="Save parameters...", uvalue="save_par")
+ menu1b99=WIDGET_BUTTON(menu1, value="Exit", uvalue="exit")
+;menu2=WIDGET_BUTTON(bar, value="Image", /MENU)
+;menu2b1=WIDGET_BUTTON(menu2, value="Options", /MENU)
+;menu2b1a=WIDGET_BUTTON(menu2b1, value="Zoom", uvalue
+
+ menu0=WIDGET_BUTTON(bar, value="Help", /MENU)
+;menu0b1=WIDGET_BUTTON(menu0, value="Help", uvalue="help")
+ menu0b2=WIDGET_BUTTON(menu0, value="Online tutorial", uvalue="web")
+ menu0b3=WIDGET_BUTTON(menu0, value="About", uvalue="about")
+
+
+ basetab=WIDGET_TAB(base_id, location=0, uvalue='')
+;--------------------------------------------------
+ baseA=WIDGET_BASE(basetab, /COL, TITLE=' INPUT ', /ALIGN_LEFT)
+;--------------------------------------------------
+ baseAA=WIDGET_BASE(baseA, /ROW, /ALIGN_LEFT)
+ base_left=WIDGET_BASE(baseAA, /COL, /ALIGN_LEFT)
+ base_img=WIDGET_BASE(base_left, /COL, /ALIGN_LEFT, /FRAME, Title="Blurred image")
+
+ base_img1=WIDGET_BASE(base_img, /ROW, /ALIGN_LEFT)
+ label = WIDGET_LABEL(base_img1, value=" Blurred image --> ")
+ button = WIDGET_BUTTON(base_img1,value="Load...", UVALUE="img_load", xsize=60, ysize=30)
+ state.id.img_info=CW_FIELD(base_img, TITLE='Information: ', $
+ VALUE=' ', /NOEDIT, /COL, YSIZE=4, xsize=50)
+ state.id.img_define = WIDGET_BUTTON(base_img, value="Define sub-domains", uvalue="img_define")
+ label= WIDGET_LABEL(base_left, value=" ")
+
+ base_psf=WIDGET_BASE(base_left, /COL, /ALIGN_LEFT, /FRAME, Title="Point Spread Functions")
+ base_psf1=WIDGET_BASE(base_psf, /ROW, /ALIGN_LEFT)
+ label = WIDGET_LABEL(base_psf1, value="Point Spread Functions --> ", /ALIGN_LEFT)
+ button = WIDGET_BUTTON(base_psf1,value="Load...", UVALUE="psf_load", xsize=60, ysize=30)
+ state.id.psf_type=CW_BGROUP(base_psf,["PSFs are stored in a single FITS file.",$
+ "PSFs are stored in separate FITS files. "], $
+ ROW=2, SET_VALUE=state.par.psf_type, UVALUE='psf_type', /EXCLUSIVE)
+ state.id.text = WIDGET_TEXT(base_psf, value=[" ", " ", " ", " "], YSIZE=4, FRAME=0)
+ state.id.psf_info=CW_FIELD(base_psf, TITLE='Information: ', $
+ VALUE=' ', /NOEDIT, /COL, YSIZE=4, xsize=50)
+ state.id.threshold = WIDGET_SLIDER(base_psf, value=fix(state.par.threshold*100), $
+ xsize=100, uvalue='threshold', $
+ title="EE Threshold (percentage)")
+
+
+ base_right=WIDGET_BASE(baseAA, /COL, /ALIGN_TOP, /FRAME)
+ state.id.ibase=WIDGET_BASE(base_right, /ALIGN_TOP)
+ a_dim=512
+ state.id.ikill=WIDGET_DRAW(state.id.ibase, XSIZE=a_dim, YSIZE=a_dim, /scroll, $
+ X_SCROLL_SIZE=a_dim, Y_SCROLL_SIZE=a_dim, KEYBOARD_EVENTS=2, $
+ /button_events, uvalue='iarea', /ALIGN_CENTER, $
+ TOOLTIP='Load an image...', RETAIN=1)
+
+ baseAB=WIDGET_BASE(baseA, /ROW, /ALIGN_LEFT, /FRAME, Title="Overlap definition")
+ base_overlap1 = WIDGET_BASE(baseAB, /COL, /ALIGN_LEFT)
+ label= WIDGET_LABEL(base_overlap1, value="Overlapping region: ")
+ label= WIDGET_LABEL(base_overlap1, value=" *) Automatically computed")
+ button = WIDGET_BUTTON(base_overlap1,value="AUTO", UVALUE='auto_overlap', xsize=70, ysize=30, /ALIGN_CENTER)
+
+ base_overlap2 = WIDGET_BASE(baseAB, /COL, /ALIGN_LEFT)
+ state.id.overlap = CW_FIELD(base_overlap2, title="*) User-defined value (pixels):", $
+ /INTEGER, uvalue='overlap', /RETURN_EVENTS, $
+ value=state.par.overlap,xsize=9,/ROW)
+ button = WIDGET_BUTTON(base_overlap2,value="CHECK", UVALUE='check_overlap', xsize=70, ysize=30, /ALIGN_CENTER)
+
+ base_overlap3 = WIDGET_BASE(baseAB, /COL, /ALIGN_LEFT)
+ label= WIDGET_LABEL(base_overlap3, value="Press REDRAW or click on the")
+ label= WIDGET_LABEL(base_overlap3, value="image to see the result.")
+ button = WIDGET_BUTTON(base_overlap3,value="REDRAW", UVALUE='redraw', xsize=70, ysize=30, /ALIGN_CENTER)
+
+
+
+;--------------------------------------------------
+ baseB=WIDGET_BASE(basetab, /COL, TITLE=' DECONVOLUTION ', /ALIGN_LEFT)
+;--------------------------------------------------
+ base_dec=WIDGET_BASE(baseB, /ROW, /ALIGN_TOP, /FRAME)
+ base_dec_left = WIDGET_BASE(base_dec, /COL, /ALIGN_LEFT, /FRAME)
+;;;;
+ base_bg = WIDGET_BASE(base_dec_left, /ROW)
+ state.id.bg_choice = CW_BGROUP(base_bg,["Costant value [ADU/pixel]", $
+ "from FITS file "], $
+ /COL, LABEL_TOP="Sky Background: ", $
+ SET_VALUE=state.par.bg_choice, UVALUE='bg_choice', /EXCLUSIVE)
+ base_bg1 = WIDGET_BASE(base_bg, /COL, /ALIGN_CENTER, /BASE_ALIGN_TOP)
+ state.id.bg=CW_FIELD(base_bg1, title="", $
+ VALUE=state.par.bg, xsize=10, /ROW, /FLOAT, /ALL_EVENTS, $
+ uvalue="bg")
+ state.id.bg_file = WIDGET_BUTTON(base_bg1,value="Load...", UVALUE="bg_file", $
+ xsize=60, ysize=30)
+ state.id.bg_info=CW_FIELD(base_bg, TITLE='Information: ', $
+ VALUE=' ', /NOEDIT, /COL, YSIZE=4, xsize=35)
+;;;;
+ state.id.sigma_ron=CW_FIELD(base_dec_left, title="Read-Out Noise value [e-/pixel]:", $
+ VALUE=state.par.sigma_ron, xsize=10, /ROW, /FLOAT, /ALL_EVENTS, $
+ uvalue="sigma_ron")
+
+ state.id.gain=CW_FIELD(base_dec_left, title="GAIN value [e-/ADU]: ", $
+ VALUE=state.par.gain, xsize=10, /ROW, /FLOAT, /ALL_EVENTS, $
+ uvalue="gain")
+;;;
+ base_dec_right = WIDGET_BASE(base_dec, /COL, /ALIGN_LEFT, /FRAME)
+ base_dec_right1 = WIDGET_BASE(base_dec_right, /ROW, /ALIGN_LEFT)
+ label=WIDGET_LABEL(base_dec_right1, value="Restoration method:")
+ method=["Richardson-Lucy (RL)", "Scaled Gradient Projection (SGP) " ]
+ state.id.method = WIDGET_COMBOBOX(base_dec_right1, value=method, uvalue='method')
+ base_dec_right2 = WIDGET_BASE(base_dec_right, /ROW, /ALIGN_LEFT)
+ label=WIDGET_LABEL(base_dec_right2, value="Stopping rule: ")
+ state.id.stop_crit = WIDGET_COMBOBOX(base_dec_right2, $
+ value=[" "], $
+ uvalue='stop_crit')
+ base_dec_right3 = WIDGET_BASE(base_dec_right, /ROW, /ALIGN_LEFT)
+ state.id.tol = CW_FIELD(base_dec_right3, title = "Tolerance: ", uvalue='tol', $
+ value=state.par.tol, xsize=15, /ROW, /FLOAT, /ALL_EVENTS)
+ base_dec_right4 = WIDGET_BASE(base_dec_right, /ROW, /ALIGN_LEFT)
+ state.id.mask_choice = CW_BGROUP(base_dec_right4, ["Apply a mask to the result of each iteration"], $
+ uvalue='mask_choice', set_value=state.par.mask_choice, $
+ /NONEXCLUSIVE)
+ state.id.mask_filename = WIDGET_BUTTON(base_dec_right4,value="Load...", UVALUE="mask_filename", $
+ xsize=60, ysize=30)
+ state.id.mask_info = CW_FIELD(base_dec_right, TITLE='Information: ', $
+ VALUE=' ', /NOEDIT, /COL, YSIZE=2, xsize=35)
+;;;
+ base_dec2=WIDGET_BASE(baseB, /COL, /ALIGN_TOP, /FRAME)
+ state.id.it_choice = CW_BGROUP(base_dec2,["Fixed value for each sub-dom. ", $
+ "Bilinear interpolation ", $
+ "User-defined [manually edit the table]"], $
+ /ROW, LABEL_TOP="Number of iterations (NoI) for each sub-domain: ", $
+ SET_VALUE=state.par.it_choice, UVALUE='it_choice', $
+ /EXCLUSIVE, /NO_RELEASE)
+ base_iter = WIDGET_BASE(base_dec2, /ROW)
+ ;base_iter1 = WIDGET_BASE(base_iter, /COL, /ALIGN_CENTER, /BASE_ALIGN_TOP)
+ dummy = widget_label(base_iter, value=" ")
+ state.id.tot_iter=CW_FIELD(base_iter, title="", $
+ VALUE=state.par.tot_iter, xsize=8, /COL, /LONG, /ALL_EVENTS, $
+ uvalue="tot_iter")
+ dummy = widget_label(base_iter, value=" Corner values: ")
+ state.id.it4=WIDGET_TABLE(base_iter, value=rotate(state.par.it4,7),COLUMN_LABELS=["Left","Right"], $
+ ROW_LABELS=["Top","Bottom"], $
+ /EDITABLE, uvalue="it4",COLUMN_WIDTH=40, ROW_HEIGHTS=27, $
+ scr_xsize=155, scr_ysize=90)
+ base_iter2 = WIDGET_BASE(base_dec2, /ROW)
+ dummy = widget_label(base_iter2, value=" ")
+ state.id.itedit=WIDGET_BUTTON(base_iter2, value="Show/Edit the NoI table", $
+ uvalue="itedit", xsize=180, ysize=30)
+
+;;;
+;base_out=WIDGET_BASE(baseB, /ROW, /ALIGN_LEFT, /FRAME)
+ base_run=WIDGET_BASE(baseB, /COL, /ALIGN_CENTER, /FRAME)
+ dummy=WIDGET_LABEL(base_run, value= " ")
+ state.id.run_deconv=WIDGET_BUTTON(base_run, value=" RUN DECONVOLUTION ", uvalue='go', $
+ xsize=250, ysize=50)
+ dummy=WIDGET_LABEL(base_run, value= " ")
+
+;--------------------------------------------------
+ baseC=WIDGET_BASE(basetab, /ROW, TITLE=' OUTPUT ', /ALIGN_LEFT)
+;--------------------------------------------------
+ base_left=WIDGET_BASE(baseC, /COL, /ALIGN_LEFT)
+ base_out=WIDGET_BASE(base_left, /COL, /ALIGN_LEFT, /FRAME)
+ label=WIDGET_LABEL(base_out, value="Reconstructed object:")
+ base_out2 = WIDGET_BASE(base_out, /ROW, /ALIGN_LEFT)
+ button = WIDGET_BUTTON(base_out2, value="Save as...", uvalue="save_out")
+ button = WIDGET_BUTTON(base_out2, value="Show", uvalue="show_out")
+ button = WIDGET_BUTTON(base_out2, value="Display FITS header", uvalue="header")
+ state.id.out_info=CW_FIELD(base_out, TITLE='Information: ', $
+ VALUE=' ', /NOEDIT, /COL, YSIZE=2, xsize=40)
+
+ base_res=WIDGET_BASE(base_left, /COL, /ALIGN_LEFT, /FRAME)
+ label=WIDGET_LABEL(base_res, value="Residual:")
+ base_res2 = WIDGET_BASE(base_res, /ROW, /ALIGN_LEFT)
+ button = WIDGET_BUTTON(base_res2, value="Save as...", uvalue="save_res")
+ button = WIDGET_BUTTON(base_res2, value="Show", uvalue="show_res")
+ button = WIDGET_BUTTON(base_res2, value="Statistics", uvalue="stat_res")
+ state.id.res_info=CW_FIELD(base_res, TITLE='Information: ', $
+ VALUE=' ', /NOEDIT, /COL, YSIZE=2, xsize=40)
+
+ base_right=WIDGET_BASE(baseC, /COL, /ALIGN_TOP, /FRAME)
+ state.id.obase=WIDGET_BASE(base_right, /ALIGN_TOP)
+ a_dim=512
+ state.id.okill=WIDGET_DRAW(state.id.obase, XSIZE=a_dim, YSIZE=a_dim, /scroll, $
+ X_SCROLL_SIZE=a_dim, Y_SCROLL_SIZE=a_dim, KEYBOARD_EVENTS=2, $
+ /button_events, uvalue='oarea', /ALIGN_CENTER, $
+ TOOLTIP='', RETAIN=1)
+
+; -----------------------------------------------------------------------
+;Standard buttons
+; -----------------------------------------------------------------------
+ base2 = WIDGET_BASE(base_id, FRAME=1, /ROW, /ALIGN_LEFT)
+ help_id = WIDGET_BUTTON(base2, VALUE='HELP',UVALUE='web')
+ cancel_id = WIDGET_BUTTON(base2, VALUE='EXIT',UVALUE='exit')
+;restore_id = WIDGET_BUTTON(base2, VALUE='RESTORE PARAMETERS',UVALUE='restore')
+;save_id = WIDGET_BUTTON(base2, VALUE='SAVE PARAMETERS',UVALUE='save')
+
+ WIDGET_CONTROL, cancel_id, /CANCEL_BUTTON
+;WIDGET_CONTROL, save_id, /DEFAULT_BUTTON
+
+ WIDGET_CONTROL, base_id, /REALIZE
+ WIDGET_CONTROL, state.id.ikill, GET_VALUE=dID
+ state.id.iarea=dID
+ WIDGET_CONTROL, base_id, SET_UVALUE=state
+ patch_set_gui, state
+ XMANAGER, 'patch_gui', base_id, GROUP_LEADER=group
+
+end
diff --git a/Patch_lib/patch_img_define.pro b/Patch_lib/patch_img_define.pro
new file mode 100644
index 0000000000000000000000000000000000000000..03b18631b1475eaed46213b1f8f14c2cf4c00d29
--- /dev/null
+++ b/Patch_lib/patch_img_define.pro
@@ -0,0 +1,144 @@
+pro patch_img_define_event, event
+
+COMMON image, img, header, N, M, a_dim, Nori, Mori
+common img_define, param, orig, N_ori, M_ori
+
+WIDGET_CONTROL, event.id, GET_UVALUE=uvalue
+WIDGET_CONTROL, event.top, GET_UVALUE=id
+case uvalue of
+'dom_type' : if (event.value EQ 1 ) then begin
+ widget_control, id.domx, SENSITIVE=0
+ widget_control, id.domy, SENSITIVE=0
+ widget_control, id.rows, /SENSITIVE
+ widget_control, id.cols, /SENSITIVE
+ endif else begin
+ widget_control, id.domx, /SENSITIVE
+ widget_control, id.domy, /SENSITIVE
+ widget_control, id.rows, SENSITIVE=0
+ widget_control, id.cols, SENSITIVE=0
+ endelse
+
+'rows' : begin
+ if (M mod event.value NE 0) then begin
+ a=dialog_message(["The chosen number of rows does not fit the image size,", $
+ "do you want to apply a zero-padding around the image?"], $
+ title="QUESTION", /QUESTION, /DEFAULT_NO, $
+ dialog_parent=event.top)
+ if (a EQ 'Yes') then begin
+ param.rows=event.value
+ new_dim = ceil(double(M) / double(param.rows)) * param.rows
+ ok=dialog_message(["Previous size = "+strtrim(N,1)+" x "+strtrim(M,1), $
+ "New size = "+strtrim(N,1)+" x "+strtrim(new_dim,1)], $
+ dialog_parent=event.top)
+ M=new_dim
+ param.domy = fix(M/param.rows)
+ widget_control, id.domy, set_value=strtrim(param.domy,1)
+ endif else begin
+ ;reset the "rows" value
+ widget_control, id.rows, set_value=strtrim(param.rows,1)
+ endelse
+ endif else begin
+ param.rows=event.value
+ param.domy = fix(M/param.rows)
+ widget_control, id.domy, set_value=strtrim(param.domy,1)
+ endelse
+ end
+'cols' : begin
+ if (N mod event.value NE 0) then begin
+ a=dialog_message(["The chosen number of columns does not fit the image size,", $
+ "do you want to apply a zero-padding around the image?"], $
+ title="QUESTION", /QUESTION, /DEFAULT_NO, $
+ dialog_parent=event.top)
+ if (a EQ 'Yes') then begin
+ param.cols=event.value
+ new_dim = ceil(double(N)/ double(param.cols)) * param.cols
+ ok=dialog_message(["Previous size = "+strtrim(N,1)+" x "+strtrim(M,1), $
+ "New size = "+strtrim(new_dim,1)+" x "+strtrim(M,1)], $
+ dialog_parent=event.top)
+ N=new_dim
+ param.domx = fix(N/param.cols)
+ widget_control, id.domx, set_value=strtrim(param.domx,1)
+ endif else begin
+ ;reset the "cols" value
+ widget_control, id.cols, set_value=strtrim(param.cols,1)
+ endelse
+ endif else begin
+ param.cols=event.value
+ param.domx = fix(N/param.cols)
+ widget_control, id.domx, set_value=strtrim(param.domx,1)
+ endelse
+ end
+;'overlap' : param.overlap=event.value
+'reset': begin
+ N=N_ori
+ M=M_ori
+ param.domx=orig.domx
+ param.domy=orig.domy
+ param.rows=orig.rows
+ param.cols=orig.cols
+ ;param.overlap=orig.overlap
+ widget_control, id.rows, set_value=param.rows
+ widget_control, id.cols, set_value=param.cols
+ widget_control, id.domx, set_value=strtrim(param.domx,1)
+ widget_control, id.domy, set_value=strtrim(param.domy,1)
+ end
+'ok' : begin
+ WIDGET_CONTROL, event.top, /DESTROY
+ end
+endcase
+
+end
+
+
+;;;;;;;;;;;;;;;;;;;;;;;;;
+; GUI generation code
+;;;;;;;;;;;;;;;;;;;;;;;;;
+function patch_img_define, title, state, GROUP_LEADER=group
+
+common image, img, header, N, M
+common img_define, param, orig, N_ori, M_ori
+
+;copy initial params in order to reset the GUI values
+param=state.par
+orig=state.par
+N_ori=N
+M_ori=M
+
+;param is the structure (of the parameters) containing edited/modified values
+;orig is the structure containing original values
+
+id = { root_base_id0: 0L, rows:0L, cols:0L, domx:0L, domy:0L, overlap:0L }
+
+modal = n_elements(group) ne 0
+root_base_id0 = widget_base(TITLE=title,/COL,MODAL=modal, GROUP_LEADER=group)
+base = widget_base(root_base_id0, /FRAME,/COL)
+label=WIDGET_LABEL(base, value="Enter the number of sub-domains (and PSFs)")
+label=WIDGET_LABEL(base, value="in the form of 2D matrix (rows x columns) ")
+label=WIDGET_LABEL(base, value=" -- press ENTER after setting each value -- ")
+base2=WIDGET_BASE(base, /ROW, /ALIGN_CENTER)
+id.rows=CW_FIELD(base2, title="", value=param.rows, /INTEGER, uvalue="rows", $
+ /RETURN_EVENTS, xsize=5 )
+label=WIDGET_LABEL(base2, value = " x ")
+id.cols=CW_FIELD(base2, title="", value=param.cols, /INTEGER, uvalue="cols", $
+ /RETURN_EVENTS, xsize=5 )
+label=WIDGET_LABEL(base, value="Resulting size of each sub-domain:")
+base3=WIDGET_BASE(base, /ROW, /ALIGN_CENTER)
+id.domx=WIDGET_TEXT(base3, value=strtrim(param.domx,1))
+label=WIDGET_LABEL(base3, value = " x ")
+id.domy=WIDGET_TEXT(base3, value=strtrim(param.domy,1))
+;label=WIDGET_LABEL(base, value=" --- ")
+;id.overlap=CW_FIELD(base, title="Overlapping region (pixels):", $
+; /INTEGER, uvalue="overlap", /RETURN_EVENTS, value=param.overlap)
+; standard buttons base
+btn_base_id = widget_base(root_base_id0,/FRAME,/ROW, /ALIGN_RIGHT)
+reset_id = widget_button(btn_base_id, VALUE="Reset to default settings", UVALUE="reset" )
+ok_id = widget_button(btn_base_id, VALUE=" Save and close ", UVALUE="ok" )
+
+
+
+WIDGET_CONTROL, root_base_id0, SET_UVALUE=id
+WIDGET_CONTROL, root_base_id0, /REALIZE
+
+xmanager, 'patch_img_define', root_base_id0, /NO_BLOCK, GROUP_LEADER=group
+return, param
+end
diff --git a/Patch_lib/patch_img_options.pro b/Patch_lib/patch_img_options.pro
new file mode 100644
index 0000000000000000000000000000000000000000..e514bdef7910a542c40609259a8e929e914167b3
--- /dev/null
+++ b/Patch_lib/patch_img_options.pro
@@ -0,0 +1,92 @@
+pro patch_img_options_event, event
+
+;COMMON image, img, header, N, M
+common img_options, param;, orig, N_ori, M_ori
+
+WIDGET_CONTROL, event.id, GET_UVALUE=uvalue
+WIDGET_CONTROL, event.top, GET_UVALUE=id
+case uvalue of
+'zoom': param.img_zoom=event.value
+'z1': begin
+ param.img_zoom=0.25
+ widget_control, id.zoom, set_value=param.img_zoom
+ end
+'z2': begin
+ param.img_zoom=0.5
+ widget_control, id.zoom, set_value=param.img_zoom
+ end
+'z3': begin
+ param.img_zoom=1.
+ widget_control, id.zoom, set_value=param.img_zoom
+ end
+'z4': begin
+ param.img_zoom=2.
+ widget_control, id.zoom, set_value=param.img_zoom
+ end
+'z5': begin
+ param.img_zoom=4.
+ widget_control, id.zoom, set_value=param.img_zoom
+ end
+'scale':param.img_scale=event.index
+'color':param.img_color=event.index
+'ok' : begin
+ WIDGET_CONTROL, event.top, /DESTROY
+ end
+endcase
+;
+end
+
+
+;;;;;;;;;;;;;;;;;;;;;;;;;
+; GUI generation code
+;;;;;;;;;;;;;;;;;;;;;;;;;
+function patch_img_options, title, state, GROUP_LEADER=group
+
+;common image, img, header, N, M
+common img_options, param;, orig, N_ori, M_ori
+
+;copy initial params in order to reset the GUI values
+param=state.par
+;orig=state.par
+;N_ori=N
+;M_ori=M
+
+;param is the structure (of the parameters) containing edited/modified values
+;orig is the structure containing original values
+
+id = { root_base_id0: 0L, zoom:0L, scale:0L, color:0L }
+
+;modal = n_elements(group) ne 0
+root_base_id0 = widget_base(TITLE=title,/COL, /modal, GROUP_LEADER=group)
+base = widget_base(root_base_id0, /FRAME,/COL)
+id.zoom=CW_FIELD(base, title="Zoom [ENTER]: ", value=param.img_zoom, /FLOATING, uvalue="zoom", $
+ /RETURN_EVENTS, xsize=10 )
+base1 =widget_base(base, /ROW)
+button1=WIDGET_BUTTON(base1, value='0.25 x', xsize=50, uvalue="z1")
+button2=WIDGET_BUTTON(base1, value='0.5 x ', xsize=50, uvalue="z2")
+button3=WIDGET_BUTTON(base1, value=' 1 x ', xsize=50, uvalue="z3")
+button4=WIDGET_BUTTON(base1, value=' 2 x ', xsize=50, uvalue="z4")
+button5=WIDGET_BUTTON(base1, value=' 4 x ', xsize=50, uvalue="z5")
+
+scale_str=['Linear', 'Log', 'Square', 'Sqrt', 'Gamma (0.2)', 'Gamma (0.1)']
+label=WIDGET_LABEL(base, value="Scale : ")
+id.scale=WIDGET_COMBOBOX(base, value=scale_str, uvalue='scale')
+
+loadct, GET_NAMES=color_str, /SILENT
+label=WIDGET_LABEL(base, value="Color : ")
+id.color=WIDGET_COMBOBOX(base, value=color_str, uvalue='color')
+
+
+; standard buttons base
+btn_base_id = widget_base(root_base_id0,/FRAME,/ROW, /ALIGN_RIGHT)
+ok_id = widget_button(btn_base_id, VALUE=" Apply and close ", UVALUE="ok")
+
+
+
+WIDGET_CONTROL, root_base_id0, SET_UVALUE=id
+WIDGET_CONTROL, root_base_id0, /REALIZE
+WIDGET_CONTROL, id.scale, set_combobox_select=param.img_scale
+WIDGET_CONTROL, id.color, set_combobox_select=param.img_color
+xmanager, 'patch_img_options', root_base_id0, /NO_BLOCK, GROUP_LEADER=group
+return, param
+end
diff --git a/Patch_lib/patch_load_background.pro b/Patch_lib/patch_load_background.pro
new file mode 100644
index 0000000000000000000000000000000000000000..6df81d48dd9d9066e50c6249a9ae32caef27ca25
--- /dev/null
+++ b/Patch_lib/patch_load_background.pro
@@ -0,0 +1,58 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+function patch_load_background, state, dummy, RESET=reset, VERBOSE=verbose
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+COMMON image, img, header, N, M, a_dim, Nori, Mori
+COMMON background, bg
+
+if KEYWORD_SET(reset) THEN reset=1 ELSE reset=0
+if KEYWORD_SET(verbose) then verb=1 else verb=0
+
+bg=readfits(state.par.bg_file, /SILENT)
+
+;if (size(bg))[0] EQ 2 then P=1
+if ((size(bg))[0] EQ 3 ) then begin
+ if (verb) then a=dialog_message(["Multiple images are not supported", $
+ "in this version of PATCH."], $
+ dialog_parent=dummy, /ERROR)
+ return, 0
+endif
+if (n_elements(bg) GT n_elements(img)) then begin
+ if (verb) then a=dialog_message(["The background array must be of the ",$
+ "same dimensions of the input image."],$
+ dialog_parent=dummy, /ERROR)
+ return, 0
+endif
+if (n_elements(bg) LT n_elements(img)) then begin
+ if (verb) then begin
+ ;ask for zero-padding
+ a=dialog_message(["The background array does not fit the image size,", $
+ "do you want to apply a zero-padding around the image?"], $
+ title="QUESTION", /QUESTION, /DEFAULT_NO, $
+ dialog_parent=dummy)
+ if (a EQ 'Yes') then bg = zero_padding(bg, N, M) else return, 0
+ endif else begin
+ ;in this case PATCH will apply zero-padding by default
+ bg = zero_padding(bg, N, M)
+ endelse
+endif
+
+if (reset) then begin
+;nothing to do!
+endif
+
+;add moment of the array
+
+
+if (verb) then begin
+ info=["filename: "+state.par.bg_file, $
+ "size: "+strtrim(N,1)+' x '+strtrim(M,1), $
+ "Mean, Variance = "+strtrim((moment(bg))[0],1)+$
+ ', '+strtrim((moment(bg))[1],1), $
+ "Skewness, Kurtosis = "+strtrim((moment(bg))[2],1)+$
+ ', '+strtrim((moment(bg))[3],1)]
+ ;show info
+ WIDGET_CONTROL, state.id.bg_info, SET_VALUE=info
+endif
+return, 1
+end
+
diff --git a/Patch_lib/patch_load_image.pro b/Patch_lib/patch_load_image.pro
new file mode 100644
index 0000000000000000000000000000000000000000..d4e584a9565bf4d6bab4febdc4e89a65a7ce27cc
--- /dev/null
+++ b/Patch_lib/patch_load_image.pro
@@ -0,0 +1,40 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+function patch_load_image, state, dummy, RESET=reset, VERBOSE=verbose
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+COMMON image, img, header, N, M, a_dim, Nori, Mori
+
+if KEYWORD_SET(reset) THEN reset=1 ELSE reset=0
+if KEYWORD_SET(verbose) then verb=1 else verb=0
+
+img=readfits(state.par.img_filename, header, /SILENT)
+
+N=(size(img))[1]
+M=(size(img))[2]
+Nori=N
+Mori=M
+if (size(img))[0] EQ 2 then P=1
+if ((size(img))[0] EQ 3 ) AND (verb) then begin
+ a=dialog_message(["Multiple images are not supported", $
+ "in this version of PATCH."], $
+ dialog_parent=dummy, /ERROR)
+ return, 0
+endif
+if (reset) then begin
+ state.par.cols=1
+ state.par.rows=1
+ state.par.overlap=0
+ state.par.domx=N/state.par.cols
+ state.par.domy=M/state.par.rows
+endif
+if (verb) then begin
+ info=["filename: "+state.par.img_filename, $
+ "size: "+strtrim(N,1)+' x '+strtrim(M,1), $
+ "domains ("+strtrim(state.par.rows*state.par.cols,1)+"): size "+$
+ strtrim(state.par.domx,1)+' x '+strtrim(state.par.domy,1), $
+ "current domain: none"]
+ ;show info
+ WIDGET_CONTROL, state.id.img_info, SET_VALUE=info
+endif
+return, 1
+end
+
diff --git a/Patch_lib/patch_load_psfs.pro b/Patch_lib/patch_load_psfs.pro
new file mode 100644
index 0000000000000000000000000000000000000000..4c0f1282f17d47a046519f443aa958641f700f50
--- /dev/null
+++ b/Patch_lib/patch_load_psfs.pro
@@ -0,0 +1,90 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+function patch_load_psfs, state, dummy
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+nb_of_psfs = state.par.cols*state.par.rows
+ ;PSF TYPE = 0 --> single file (cube)
+ if state.par.psf_type EQ 0 then begin
+ psf = readfits((*state.par.psf_filename)[0], /SILENT)
+ NP=(size(psf))[1]
+ MP=(size(psf))[2]
+ if (size(psf))[0] EQ 2 then PP=1 else PP=(size(psf))[3]
+ if (PP NE nb_of_psfs) then begin
+ a=dialog_message(["You must load a cube FITS file. The ", $
+ "dimensions should be NxMxP, where P ",$
+ "is the number of sub-domains."], $
+ dialog_parent=dummy, /ERROR)
+ return, 0
+ endif
+ if (PP EQ 1) then begin
+ info=["filename: "+(*state.par.psf_filename)[0], $
+ "size: "+strtrim(NP,1)+' x '+strtrim(MP,1), $
+ " ", $
+ " "]
+ endif else begin
+ info=["filename: "+(*state.par.psf_filename)[0], $
+ "size: "+strtrim(NP,1)+' x '+strtrim(MP,1)+' x '+strtrim(PP,1), $
+ " ", $
+ " "]
+ endelse
+ ;show info
+ WIDGET_CONTROL, state.id.psf_info, SET_VALUE=info
+ if (NP LT state.par.domx+2*state.par.overlap) OR $
+ (MP LT state.par.domy+2*state.par.overlap) then begin
+ a=dialog_message(["The dimensions of the PSFs will be zero-padded", $
+ "in order to fit the size of the sub-domains. ",$
+ "(plus the 2*overlap value)."], $
+ dialog_parent=dummy)
+ endif
+ endif
+ ; PSF TYPE = 1 --> multiple files
+ if state.par.psf_type EQ 1 then begin
+ psf0=readfits((*state.par.psf_filename)[0], /SILENT)
+ NP=(size(psf0))[1]
+ MP=(size(psf0))[2]
+ if (size(psf0))[0] EQ 2 then PP=1 else begin
+ a=dialog_message("Please, check the dimensions of the first PSF. ", $
+ dialog_parent=dummy, /ERROR)
+ return, 0
+ endelse
+ info=["filename: "+(*state.par.psf_filename)[0], $
+ "size: "+strtrim(NP,1)+' x '+strtrim(MP,1), $
+ "First PSF loaded... ", $
+ "... checking for other files"]
+ ;show info (first PSF)
+ WIDGET_CONTROL, state.id.psf_info, SET_VALUE=info
+ ok_go_on=1
+ for k=1,nb_of_psfs-1 do begin
+ psf=readfits((*state.par.psf_filename)[k],/SILENT)
+ if (size(psf0, /N_ELE) NE size(psf, /N_ELE)) then begin
+ a=dialog_message("Please, check the dimensions of PSF nb "+strtrim(k,1), $
+ dialog_parent=dummy, /ERROR)
+ ok_go_on=0
+ endif
+ endfor
+ if (NP LT state.par.domx+2*state.par.overlap) OR $
+ (MP LT state.par.domy+2*state.par.overlap) then begin
+ a=dialog_message(["The dimensions of the PSFs will be zero-padded", $
+ "in order to fit the size of the sub-domains. ",$
+ "(plus the 2*overlap value)."], $
+ dialog_parent=dummy)
+ endif
+ if (ok_go_on) then begin
+ info=["filename: "+(*state.par.psf_filename)[0], $
+ "size: "+strtrim(NP,1)+' x '+$
+ strtrim(MP,1)+' x '+$
+ strtrim(nb_of_psfs,1), $
+ " ", $
+ "All files are OK."]
+ endif else begin
+ info=["filename: "+(*state.par.psf_filename)[0], $
+ "size: "+strtrim(NP,1)+' x '+$
+ strtrim(MP,1), $
+ " ERROR! ", $
+ "Check your files and re-load PSFs "]
+ endelse
+ ;show info (OK or ERROR)
+ WIDGET_CONTROL, state.id.psf_info, SET_VALUE=info
+ return, 1
+ endif
+end
diff --git a/Patch_lib/patch_noi_table.pro b/Patch_lib/patch_noi_table.pro
new file mode 100644
index 0000000000000000000000000000000000000000..2b59237480d20ceb1021f630b89a05e514921632
--- /dev/null
+++ b/Patch_lib/patch_noi_table.pro
@@ -0,0 +1,82 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+pro patch_noi_table_event, event
+ common table, param, orig
+
+ WIDGET_CONTROL, event.id, GET_UVALUE=uvalue
+ WIDGET_CONTROL, event.top, GET_UVALUE=id
+ case uvalue of
+ 'edit': begin
+ widget_control, id.table, get_value=dummy
+ if dummy[event.x,event.y] LE 0 then begin
+ a=dialog_message("Insert positive (>0) values!", /ERROR)
+ endif
+ dummy=rotate(dummy,7)
+ *param.tab_iter=dummy
+ end
+ 'reset': begin
+ (*param.tab_iter)[*]=orig[*]
+ dummy=orig
+ dummy=reform(dummy,param.cols,param.rows)
+ dummy=rotate(dummy,7)
+ widget_control, id.table, set_value=dummy
+ end
+ 'ok' : begin
+ WIDGET_CONTROL, event.top, /DESTROY
+ end
+
+ endcase
+
+end
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+function patch_noi_table, par, GROUP_LEADER=group
+
+ common table, param, orig
+ param=par
+ orig=(*par.tab_iter)[*]
+;param is the structure (of the parameters) containing edited/modified values
+;orig is the structure containing original values
+
+ id = { $
+ root_base_id0 : 0L, $
+ table:0L $
+ }
+ param=par
+
+ modal = n_elements(group) ne 0
+ root_base_id0 = widget_base(TITLE="Number of Iterations (Noi) Table",/COL, $
+ modal=modal, GROUP_LEADER=group)
+ base = widget_base(root_base_id0, /FRAME,/COL)
+
+ tab_val=make_array(par.cols,par.rows,/LONG)
+ ;help, tab_val
+ for k=0, par.cols*par.rows-1 do begin
+ tab_val[k mod par.cols,k/par.cols]=(*param.tab_iter)[k]
+ endfor
+ tab_val=rotate(tab_val,7)
+ id.table= WIDGET_TABLE(base, value=tab_val,$
+ ; COLUMN_LABELS=["Left","Right"], $
+ ; ROW_LABELS=["Top","Bottom"], $
+ /EDITABLE, uvalue="edit", /NO_HEADERS, $
+ COLUMN_WIDTH=40, ROW_HEIGHTS=27, $
+ scr_xsize=400, scr_ysize=300)
+
+
+
+; standard buttons base
+ btn_base_id = widget_base(root_base_id0,/FRAME,/ROW, /ALIGN_RIGHT)
+ reset_id = widget_button(btn_base_id, $
+ VALUE="Reset to default settings", $
+ UVALUE="reset" )
+
+ ok_id = widget_button(btn_base_id, $
+ VALUE=" Save and close ", $
+ UVALUE="ok" )
+
+
+
+ WIDGET_CONTROL, root_base_id0, SET_UVALUE=id
+ WIDGET_CONTROL, root_base_id0, /REALIZE
+ xmanager, 'patch_noi_table', root_base_id0, /NO_BLOCK, GROUP_LEADER=group
+ return, param
+
+end
diff --git a/Patch_lib/patch_noi_table_compute.pro b/Patch_lib/patch_noi_table_compute.pro
new file mode 100644
index 0000000000000000000000000000000000000000..9126ca05d43f27bebdaa7ad5a34962a65a8b5f84
--- /dev/null
+++ b/Patch_lib/patch_noi_table_compute.pro
@@ -0,0 +1,30 @@
+pro patch_noi_table_compute, choice, par
+
+ case choice of
+ 0: begin
+ *par.tab_iter=make_array(par.cols*par.rows, value=LONG(par.tot_iter))
+ end
+ 1: begin
+ N=par.cols*par.rows
+ if N GT 1 then begin
+ ;interpolation
+ a=double(par.it4)
+ *par.tab_iter=make_array(N,/LONG)
+ for k=0,N-1 do begin
+ i=(k mod par.cols)
+ j=(k / par.cols)
+ d1=a[0,0]+(a[1,0]-a[0,0])/(par.cols-1)*i
+ d2=a[0,1]+(a[1,1]-a[0,1])/(par.cols-1)*i
+ d3=round(d1+(d2-d1)/(par.rows-1)*j)
+ (*par.tab_iter)[k]=d3
+ endfor
+ endif else (*par.tab_iter)[*]=1L
+ end
+ 2: begin
+ if n_elements(*par.tab_iter) NE par.cols*par.rows then begin
+ *par.tab_iter=make_array(par.cols*par.rows, value=1L)
+ endif
+ end
+ endcase
+ ;print, "--> computed: ", *par.tab_iter
+end
diff --git a/Patch_lib/patch_psf_ee.pro b/Patch_lib/patch_psf_ee.pro
new file mode 100644
index 0000000000000000000000000000000000000000..f147e7dce5c764c4a02fc39221f034e751dc61dd
--- /dev/null
+++ b/Patch_lib/patch_psf_ee.pro
@@ -0,0 +1,40 @@
+;Last update: 2014-10-16 by Andrea La Camera
+;+
+;-
+
+function patch_psf_ee, psf, threshold
+on_error,2
+
+t0=systime(1)
+if (threshold LT 0.) or (threshold GT 1.) then message, 'Threshold must be in range [0,1]'
+
+dim=max(size(psf, /DIM))
+
+if (dim mod 2) EQ 1 then dim++
+
+psf2=zero_padding(psf,dim,dim)
+ee=dblarr(dim/2) ;Enclosed Energy as a function of the radius (distance from
+ ;the centre of the PSF)
+stpcrt=1
+k0=0.
+k1=dim/2.
+k=0
+if threshold GT 0 then begin
+ while (stpcrt) do begin
+ k=ceil((k0+k1)/2.)
+ ee[k-1]=total(psf[dim/2-k:dim/2+k-1, dim/2-k:dim/2+k-1],/DOUBLE)
+ if ee[k-1] GT threshold then begin
+ k1=k
+ endif else begin
+ k0=k
+ endelse
+ if abs(k0-k1) LE 1 then stpcrt = 0
+ endwhile
+endif
+;plot, ee
+;oplot, make_array(dim/2, value=threshold)
+;print, 'Execution time = ', systime(1)-t0
+;print, "Delta value = ", k
+return, k
+
+end
diff --git a/Patch_lib/patch_set_gui.pro b/Patch_lib/patch_set_gui.pro
new file mode 100644
index 0000000000000000000000000000000000000000..6fbef9691454d6774d39cd77b7ff6c294fcaa340
--- /dev/null
+++ b/Patch_lib/patch_set_gui.pro
@@ -0,0 +1,80 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+pro patch_set_gui, state
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+if state.par.psf_type EQ 0 then begin
+ value=["The file dimensions must be NxMxP, where P is the",$
+ "number of sub-domains. The PSFs must be ordered ",$
+ "from 0 to P-1, as shown in the image. ", $
+ " "]
+endif else begin
+ value=["The files (one for each sub-domain) must have a ",$
+ "common prefix, followed by a number indicating the",$
+ "corresponding sub-domain. E.g. (16 sub-domains): ",$
+ "'psf_00.fits', 'psf_01.fits', ..., 'psf_15.fits' "]
+endelse
+widget_control, state.id.text, set_value=value
+WIDGET_CONTROL, state.id.method, SET_COMBOBOX_SELECT=state.par.method
+
+
+if state.par.bg_choice EQ 0 then begin
+ widget_control, state.id.bg, /SENSITIVE
+ widget_control, state.id.bg_file, SENSITIVE=0
+ widget_control, state.id.bg_info, SET_VALUE=strarr(4)
+endif else begin
+ widget_control, state.id.bg, SENSITIVE=0
+ widget_control, state.id.bg_file, /SENSITIVE
+endelse
+
+
+case state.par.method of
+0: begin ;RL parameters
+ widget_control, state.id.stop_crit, SET_VALUE=['iter > MAX ITER']
+ widget_control, state.id.stop_crit, SET_COMBOBOX_SELECT=state.par.stop_crit<0
+ ;force stop_crit = 0
+ state.par.stop_crit = 0
+ end
+1: begin ;SGP parameters
+ widget_control, state.id.stop_crit, SET_VALUE=['iter > MAX ITER', $
+ 'distance of successive iterations', $
+ 'data-fidelity function', $
+ 'discrepancy principle']
+ widget_control, state.id.stop_crit, SET_COMBOBOX_SELECT=state.par.stop_crit<3
+ end
+endcase
+
+if state.par.stop_crit GT 0 then widget_control, state.id.tol, /SENSITIVE else $
+ widget_control, state.id.tol, SENSITIVE=0
+
+case state.par.mask_choice of
+ 0: begin
+ widget_control, state.id.mask_filename, SENSITIVE=0
+ widget_control, state.id.mask_info, SENSITIVE=0
+ end
+ 1: begin
+ widget_control, state.id.mask_filename, /SENSITIVE
+ widget_control, state.id.mask_info, /SENSITIVE
+ end
+endcase
+
+case state.par.it_choice of
+ 0: begin
+ widget_control, state.id.tot_iter, /SENSITIVE
+ widget_control, state.id.it4, SENSITIVE=0
+ end
+ 1: begin
+ widget_control, state.id.tot_iter, SENSITIVE=0
+ widget_control, state.id.it4, /SENSITIVE
+ if state.par.cols*state.par.rows LT 4 then begin
+ a=dialog_message("At least 4 sub-domains are required!", /ERROR)
+ endif
+ end
+ 2: begin
+ widget_control, state.id.tot_iter, SENSITIVE=0
+ widget_control, state.id.it4, SENSITIVE=0
+ end
+endcase
+
+
+
+end
diff --git a/Patch_lib/patch_show_image.pro b/Patch_lib/patch_show_image.pro
new file mode 100644
index 0000000000000000000000000000000000000000..7a1947eff79894bf10b92dd30c9b9a30ecd31193
--- /dev/null
+++ b/Patch_lib/patch_show_image.pro
@@ -0,0 +1,71 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+pro patch_show_image, state, REDRAW=redraw
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+COMMON image, img, header, N, M, a_dim, Nori, Mori
+
+NN=N+2*state.par.overlap
+MM=M+2*state.par.overlap
+
+loadct, state.par.img_color, /SILENT
+
+;show image (REDRAW if requested)
+IF KEYWORD_SET(redraw) THEN BEGIN
+ Widget_Control, state.id.ibase, Map=0
+ Widget_Control, state.id.ikill, /destroy
+ if state.par.cols*state.par.rows GT 1 then begin
+ tooltip_str='left click to select, right click for options'
+ endif else tooltip_str='right click for options'
+ drawID = Widget_Draw(state.id.ibase, /scroll, KEYBOARD_EVENTS=2, $
+ /button_events, uvalue='iarea', /ALIGN_CENTER, $
+ XSize=NN*state.par.img_zoom, $
+ YSize=MM*state.par.img_zoom, $
+ X_Scroll_Size= a_dim,$
+ Y_Scroll_Size= a_dim,$
+ TOOLTIP=tooltip_str)
+ Widget_Control, drawID, /Realize
+ Widget_Control, drawID, Get_Value=wid
+ state.id.iarea = wid
+ state.id.ikill = drawID
+ Widget_Control, state.id.ibase, Map=1
+ENDIF
+
+wset, state.id.iarea
+
+case state.par.img_scale of
+ '0': tvscl, CONGRID(zero_padding(img,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)
+ '1': begin
+ img2=LOGSCL(img, EXPONENT=.5) ;(1./(1. + (0.5/(img>1.e-5))^2.)) ;img+(1e-5*MAX(img))+ABS(min(img))
+ tvscl,CONGRID(zero_padding(img2,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)
+ end
+ '2': tvscl, DOUBLE(CONGRID(zero_padding(img,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)^2.)
+ '3': tvscl, SQRT(CONGRID(zero_padding(img,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom))
+ '4': tvscl, GMASCL(CONGRID(zero_padding(img,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom),$
+ gamma=0.2)
+ '5': tvscl, GMASCL(CONGRID(zero_padding(img,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom),$
+ gamma=0.1)
+
+endcase
+if state.par.cols*state.par.rows GT 1 then begin
+ for x=0,state.par.cols do plots, (x*state.par.domx-1+state.par.overlap)*state.par.img_zoom, $
+ (indgen(M)+state.par.overlap)*state.par.img_zoom, $
+ /DEVICE, color='FFFFFF'x, LINESTYLE=2, thick=2.0
+ for y=0,state.par.rows do plots, (indgen(N)+state.par.overlap)*state.par.img_zoom, $
+ (y*state.par.domy-1+state.par.overlap)*state.par.img_zoom, $
+ /DEVICE, color='FFFFFF'x, LINESTYLE=2, thick=2.0
+
+ if ((state.par.domx+state.par.domy)*state.par.img_zoom) GE 51 then begin
+ ; for k=0,state.par.cols*state.par.rows-1 do begin
+ ; xyouts, ((((k mod state.par.cols)+0.5)*state.par.domx)+state.par.overlap)*state.par.img_zoom, $
+ ; ((((k / state.par.cols)+0.5)*state.par.domy)+state.par.overlap)*state.par.img_zoom, $
+ ; strtrim(k,1), /DEVICE, charsize=1.5, ALIGNMENT=0.5, charthick=1.
+ ; endfor
+ arr=indgen(state.par.cols*state.par.rows)
+ xyouts, ((((arr mod state.par.cols)+0.5)*state.par.domx)+state.par.overlap)*state.par.img_zoom, $
+ ((((arr / state.par.cols)+0.47)*state.par.domy)+state.par.overlap)*state.par.img_zoom, $
+ strtrim(indgen(state.par.cols*state.par.rows),1), /DEVICE, charsize=1.5, ALIGNMENT=0.5, charthick=1.
+ endif
+endif
+
+return
+end
+
diff --git a/Patch_lib/patch_show_out.pro b/Patch_lib/patch_show_out.pro
new file mode 100644
index 0000000000000000000000000000000000000000..e484dc313fdc1a09937a3efe1f73931af33ba7ee
--- /dev/null
+++ b/Patch_lib/patch_show_out.pro
@@ -0,0 +1,52 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+pro patch_show_out, state, REDRAW=redraw
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+COMMON image, img, header, N, M, a_dim, Nori, Mori
+COMMON output, out, res_array, residual
+
+NN=Nori;N+2*state.par.overlap
+MM=Mori;M+2*state.par.overlap
+
+loadct, state.par.img_color, /SILENT
+
+;show image (REDRAW if requested)
+IF KEYWORD_SET(redraw) THEN BEGIN
+ Widget_Control, state.id.obase, Map=0
+ Widget_Control, state.id.okill, /destroy
+ if state.par.cols*state.par.rows GT 1 then tooltip_str='right click for options' $
+ else tooltip_str='right click for options'
+ drawID = Widget_Draw(state.id.obase, /scroll, KEYBOARD_EVENTS=2, $
+ /button_events, uvalue='oarea', /ALIGN_CENTER, $
+ XSize=NN*state.par.img_zoom, $
+ YSize=MM*state.par.img_zoom, $
+ X_Scroll_Size= a_dim,$
+ Y_Scroll_Size= a_dim,$
+ TOOLTIP=tooltip_str)
+ Widget_Control, drawID, /Realize
+ Widget_Control, drawID, Get_Value=wid
+ state.id.oarea = wid
+ state.id.okill = drawID
+ Widget_Control, state.id.obase, Map=1
+ENDIF
+
+wset, state.id.oarea
+
+case state.par.img_scale of
+ '0': tvscl, CONGRID(zero_padding(out,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)
+ '1': begin
+ img2=LOGSCL(out>0., EXPONENT=.5) ;(1./(1. + (0.5/(img>1.e-5))^2.)) ;img+(1e-5*MAX(img))+ABS(min(img))
+ tvscl,CONGRID(zero_padding(img2,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)
+ end
+ '2': tvscl, DOUBLE(CONGRID(zero_padding(out,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)^2.)
+ '3': tvscl, SQRT(CONGRID(zero_padding(out,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom))
+ '4': tvscl, GMASCL(CONGRID(zero_padding(out,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom),$
+ gamma=0.2)
+ '5': tvscl, GMASCL(CONGRID(zero_padding(out,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom),$
+ gamma=0.1)
+
+endcase
+
+
+return
+end
+
diff --git a/Patch_lib/patch_show_res.pro b/Patch_lib/patch_show_res.pro
new file mode 100644
index 0000000000000000000000000000000000000000..fdaf34fb736578fec2bbbf14c10bb27e1af6822a
--- /dev/null
+++ b/Patch_lib/patch_show_res.pro
@@ -0,0 +1,50 @@
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+pro patch_show_res, state, REDRAW=redraw
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+COMMON image, img, header, N, M, a_dim, Nori, Mori
+COMMON output, out, res_array, residual
+
+NN=Nori;N+2*state.par.overlap
+MM=Mori;M+2*state.par.overlap
+
+loadct, state.par.img_color, /SILENT
+
+;show image (REDRAW if requested)
+IF KEYWORD_SET(redraw) THEN BEGIN
+ Widget_Control, state.id.obase, Map=0
+ Widget_Control, state.id.okill, /destroy
+ if state.par.cols*state.par.rows GT 1 then tooltip_str='right click for options' $
+ else tooltip_str='right click for options'
+ drawID = Widget_Draw(state.id.obase, /scroll, KEYBOARD_EVENTS=2, $
+ /button_events, uvalue='oarea', /ALIGN_CENTER, $
+ XSize=NN*state.par.img_zoom, $
+ YSize=MM*state.par.img_zoom, $
+ X_Scroll_Size= a_dim,$
+ Y_Scroll_Size= a_dim,$
+ TOOLTIP=tooltip_str)
+ Widget_Control, drawID, /Realize
+ Widget_Control, drawID, Get_Value=wid
+ state.id.oarea = wid
+ state.id.okill = drawID
+ Widget_Control, state.id.obase, Map=1
+ENDIF
+
+wset, state.id.oarea
+
+case state.par.img_scale of
+ '0': tvscl, CONGRID(zero_padding(residual,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)
+ '1': begin
+ img2=LOGSCL(residual>0., EXP=0.5) ;(1./(1. + (0.5/(img>1.e-5))^2.)) ;img+(1e-5*MAX(img))+ABS(min(img))
+ tvscl,CONGRID(zero_padding(img2,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)
+ end
+ '2': tvscl, DOUBLE(CONGRID(zero_padding(residual,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom)^2.)
+ '3': tvscl, SQRT(CONGRID(zero_padding(residual,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom))
+ '4': tvscl, GMASCL(CONGRID(zero_padding(residual,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom),$
+ gamma=0.2)
+ '5': tvscl, GMASCL(CONGRID(zero_padding(residual,NN,MM),NN*state.par.img_zoom, MM*state.par.img_zoom),$
+ gamma=0.1)
+
+endcase
+return
+end
+
diff --git a/Patch_lib/sgp_deblurring.pro b/Patch_lib/sgp_deblurring.pro
new file mode 100644
index 0000000000000000000000000000000000000000..02396baac0f7339b9ee7129b651891d59979986a
--- /dev/null
+++ b/Patch_lib/sgp_deblurring.pro
@@ -0,0 +1,589 @@
+;MODIFCA PER LA MASCHERA DEL PROIETTORE
+;pr_mask
+;Andrea La Camera, 2014/03/19
+
+;MODIFICA PER CALCOLARE IL RESIDUO E LA SUA STATISTICA
+;res --> residuo calcolato come (g-(Af+b))/sqrt(Af+b)
+;res_stat --> statistica del residuo, output di histogauss (astro-lib)
+;Andrea La Camera, 2013/08/09
+
+;MODIFICA PER GESTIRE DECONVOLUZIONE DI SOLO BACKGROUND
+;1) x=(total(gn-bn) > 0.) /n_elements(gn) * blablabla
+;2) flux = (total(total(gn,1),1) - total(total(bn,1),1)) > var.eps
+;Andrea La Camera, 2013/03/19
+
+;AGGIORNAMENTO per BACKGROUND
+;bg passa da una constante a un array delle stesse dimensioni di gn
+;Andrea La Camera, 2013/02/21
+;
+;QUESTA VERSIONE E' STATA CORRETTA IL GIORNO 19 gennaio 2012
+;da Andrea La Camera, in collaborazione con M.Prato
+;
+
+
+PRO sgp_deblurring, A, gn, x, iter, err, Discr, times, pr_mask, res, res_stat, $
+ bg = bg, bound = bound, $
+ initialization = initialization, maxit = maxit, obj = obj, $
+ stopcrit = stopcrit, tol = tol, verb = verb
+
+; SGP_deblurring - SGP algorithm for non-regularized deblurring
+; This function solves an image deblurring problem by applying the SGP
+; algorithm to the minimization of the generalized Kullback-Leibler
+; divergence with no regularization [1]:
+;
+; min KL(A*x + bg, gn)
+; x in OMEGA
+;
+; where KL(u,v) is the generalized Kullback-Leibler divergence between
+; vectors u and v, bg is the background, gn are the observed data and
+; the feasible set OMEGA is x(i) >= 0;
+;
+; [1] S. Bonettini, R. Zanella, L. Zanni,
+; "A scaled gradient projection method for constrained image deblurring",
+; Inverse Problems 25(1), 2009, January, 015002.
+;
+; SYNOPSIS
+; SGP_deblurring, A, gn, x, iter, err, Discr, times [, opts]
+;
+; MANDATORY INPUT
+; A (double array) - measuring matrix used to apply the blurring
+; operator, that is to compute A*x
+; REMARK: all comlumns of A must sum-up to 1
+; gn (double array) - measured image
+;
+; OPTIONAL INPUT
+; The following options must be provided as keyword/value pairs.
+; 'OBJ' - Exact solution, for error calculation (double array)
+; 'BG' - Background value (double array)
+; DEFAULT = dblarr(size(gn))
+; 'INITIALIZATION' - Choice for starting point:
+; 0 - all zero starting point
+; 1 - initialization with gn
+; 2 - initialization with
+; ones(size(gn))*sum(gn(:) - bg(:)) / numel(gn)
+; x0 - user-provided starting point (double array)
+; DEFAULT = 2
+; 'MAXIT' - Maximum number of iterations (integer)
+; DEFAULT = 1000
+; 'VERB' - Verbosity level (integer)
+; 0 - silent
+; 1 - print configuration parameters at startup and
+; some information at each iteration
+; DEFAULT = 0
+; 'STOPCRIT' - Choice for stopping rule (integer)
+; 1 -> iter > MAXIT
+; 2 -> ||x_k - x_(k-1)|| <= tol*||x_k|| OR iter > MAXIT
+; 3 -> |KL_k - KL_(k-1)| <= tol*|KL_k| OR iter > MAXIT
+; 4 -> (2/N)*KL_k <= tol OR iter > MAXIT
+; DEFAULT = 1;
+; 'TOL' - Tolerance used in the stopping criterion
+; DEFAULT = 1e-4 if STOPCRITERION = 2 or 3
+; DEFAULT = 1+1/mean(gn) if STOPCRITERION = 4
+; 'BOUND' - Flag to enable bound effect reduction
+;
+; OUTPUT
+; x - Reconstructed data
+; iter - Number of iterations
+; err - Error value at each iteration. If OBJ was not given,
+; then err is the empty matrix.
+; discr - Discrepancy value after each iteration:
+; D = 2/numel(x_k) * KL( Ax_k + bg, gn)
+; times - Time elapsed after each iteration
+;
+; ------------------------------------------------------------------------------
+;
+; This software is developed within the research project
+;
+; PRISMA - Optimization methods and software for inverse problems
+; http://www.unife.it/prisma
+;
+; funded by the Italian Ministry for University and Research (MIUR), under
+; the PRIN2008 initiative, grant n. 2008T5KA4L, 2010-2012. This software is
+; part of the package "IRMA - Image Reconstruction in Microscopy and Astronomy"
+; currently under development within the PRISMA project.
+;
+; Version: 1.0
+; Date: July 2011
+
+; Authors:
+; Roberto Cavicchioli, Marco Prato, Luca Zanni
+; Dept. of Pure Appl. Math., Univ. of Modena and Reggio Emilia, Italy
+; roberto.cavicchioli@unimore.it, marco.prato@unimore.it, luca.zanni@unimore.it
+; Mario Bertero, Patrizia Boccacci
+; DISI (Dipartimento di Informatica e Scienze dell'Informazione), University of Genova, Italy
+; bertero@disi.unige.it, boccacci@disi.unige.it
+;
+; Software homepage: http://www.unife.it/prisma/software
+;
+; Copyright (C) 2011 by M. Bertero, P. Boccacci, R. Cavicchioli, M. Prato, L. Zanni
+; ------------------------------------------------------------------------------
+; COPYRIGHT NOTIFICATION
+;
+; Permission to copy and modify this software and its documentation for
+; internal research use is granted, provided that this notice is retained
+; thereon and on all copies or modifications. The authors and their
+; respective Universities makes no representations as to the suitability
+; and operability of this software for any purpose. It is provided "as is"
+; without express or implied warranty. Use of this software for commercial
+; purposes is expressly prohibited without contacting the authors.
+;
+; This program is free software; you can redistribute it and/or modify it
+; under the terms of the GNU General Public License as published by the
+; Free Software Foundation; either version 3 of the License, or (at your
+; option) any later version.
+;
+; This program is distributed in the hope that it will be useful, but
+; WITHOUT ANY WARRANTY; without even the implied warranty of
+; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
+; See the GNU General Public License for more details.
+;
+; You should have received a copy of the GNU General Public License along
+; with this program; if not, either visite http://www.gnu.org/licenses/
+; or write to
+; Free Software Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
+; ==============================================================================
+
+; start the clock
+t0 = systime(2)
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;
+; SGP default parameters ;
+;;;;;;;;;;;;;;;;;;;;;;;;;;
+alpha_min = 1.e-5 ; alpha lower bound
+alpha_max = 1.e5 ; alpha upper bound
+theta = 0.4 ; backtracking parameter
+beta = 1.e-4 ; for sufficient decrease
+initalpha = 1.3 ; initial alpha
+M = 1 ; memory in obj. function value (if M = 1 monotone)
+Malpha = 3 ; alfaBB1 memory
+tau = 0.5 ; alternating parameter
+initflag = 2 ; 2 -> init with constant
+errflag = 0B ; 0 -> no error calculation
+err = -1
+
+if (keyword_set(bound) - 1) then bound = 0B ; bound effects
+if (keyword_set(maxit) - 1) then maxit = 1000 ; maximum number of iterations
+if (keyword_set(stopcrit) - 1) then stopcrit = 1 ; 1 -> number of iterations
+if (keyword_set(verb) - 1) then verb = 0 ; 0 -> silent
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+dim = size(gn,/dimensions)
+if (size(gn,/n_dimensions) eq 2) then dim = [dim,1]
+
+if (keyword_set(bg) - 1) then bg = dblarr(dim) ; background value
+
+if bound then begin
+ newdim = [2^(floor(alog(dim[0])/alog(2))+1) , 2^(floor(alog(dim[1])/alog(2))+1)]
+ margine = (newdim - dim[0:1])/2
+ ;; WARNING: at this level, if A and gn differ in dimensions, then
+ ;; the input A MUST be sized as newdim[0] x newdim[1] by the user
+ dim_A=size(A,/dimensions)
+ if n_elements(dim_A) eq 2 then dim_A=[dim_A,1]
+ if array_equal(dim_A,dim) then begin
+ psf_work = dblarr(newdim[0],newdim[1],dim[2])
+ psf_work[margine[0]:margine[0]+dim[0]-1,margine[1]:margine[1]+dim[1]-1,*] = A
+ endif else begin
+ psf_work = A
+ end
+ TF = dcomplexarr(newdim[0],newdim[1],dim[2])
+ CTF = TF
+ for i = 0, dim[2]-1 do begin
+ TF[*,*,i] = fft(shift(reform(psf_work[*,*,i]),newdim/2))*newdim[0]*newdim[1]
+ CTF[*,*,i] =conj(TF[*,*,i])
+ endfor
+endif else begin
+ TF = dcomplexarr(dim)
+ CTF = TF
+ for i = 0, dim[2]-1 do begin
+ TF[*,*,i] = fft(shift(reform(A[*,*,i]),dim[0:1]/2))*dim[0]*dim[1]
+ CTF[*,*,i] = conj(TF[*,*,i])
+ endfor
+end
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; Read the optional parameters ;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+if arg_present(obj) then errflag = 1B
+if keyword_set(initialization) then begin
+ if product(size(initialization,/dimensions)) gt 1 then begin
+ initflag = 999
+ x = initialization ; initial x provided by user
+ endif else begin
+ initflag = initialization
+ endelse
+endif
+
+;;;;;;;;;;;;;;;;;;
+; starting point ;
+;;;;;;;;;;;;;;;;;;
+case initflag of
+ 0: x = dblarr(dim[0],dim[1]) ; all zeros
+ 1: x = gn[*,*,0] ; gn
+ 2: x = (total(gn-bg)>0.)/n_elements(gn)*make_array(dim[0],dim[1],value=1.,/double) ; same flux as gn - bg
+ 999: if not(array_equal(size(x,/dimensions),dim[0:1])) then begin ; x is explicitly given, check dimensions
+ message, 'Invalid size of the initial point.'
+ endif
+ else: message, 'Unknown initialization option.'
+endcase
+; size of the images
+obj_size = dim[0:1]
+if bound then begin
+ work_size = newdim
+endif else begin
+ work_size = obj_size
+endelse
+; mask of the projector (adapting size)
+pr_mask=zero_padding(pr_mask,work_size[0], work_size[1], value=1.)
+;help, pr_mask
+;;;;;;;;;;;;;;;;;;
+; stop criterion ;
+;;;;;;;;;;;;;;;;;;
+
+if (stopcrit ne 1 and stopcrit ne 2 and stopcrit ne 3 and stopcrit ne 4) then begin
+ message, 'Unknown stopping criterion:', stopcrit
+end
+
+if (keyword_set(tol) - 1) then begin
+ if (stopcrit eq 2 or stopcrit eq 3) then tol = 1.e-4
+ if (stopcrit eq 4) then tol = 1. + 1./mean(gn)
+end
+
+;;;;;;;;;;;;;;;;
+; data scaling ;
+;;;;;;;;;;;;;;;;
+scaling = max(gn)
+gn = gn/scaling
+bg = bg/scaling
+x = x/scaling
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; change the null pixels of the observed image ;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+par = machar(/double)
+vmin = min(gn[where(gn gt 0)])
+gn_nonpos = where(gn le 0, n_nonpos)
+if (n_nonpos gt 0) then gn[gn_nonpos] = vmin*par.eps^2
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; initializations and computations that need only once ;
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+N = n_elements(gn)/dim[2] ; pixels in the image
+flux = (total(total(gn,1),1)-total(total(bg,1),1)) > par.eps ; exact flux
+iter = 1 ; iteration counter
+alpha = initalpha ; initial alpha
+Valpha = alpha_max * (dblarr(Malpha)+1.) ; memory buffer for alpha
+Fold = -1.e30 * (dblarr(M)+1.) ; memory buffer for obj. func.
+Discr_coeff = 2./(N*dim[2])*scaling ; discrepancy coefficient
+ONE = make_array(work_size,value = dim[2], /double)
+
+if bound then begin
+ Mask = dblarr(work_size)
+ Mask[margine[0]:margine[0]+dim[0]-1,margine[1]:margine[1]+dim[1]-1] = 1.
+ work_index = long(where(Mask gt 0.))
+ sigma = 1.e-2
+ Weight = afunction(Mask,CTF)
+ no_z = where(Weight[*,*,0] gt sigma)
+ if dim[2] gt 1 then begin
+ tmp = intarr(newdim[0],newdim[1])
+ tmp(no_z)=1
+ for i = 1, dim[2]-1 do begin
+ no_z_tmp = where(Weight[*,*,i] gt sigma)
+ tmp2 = intarr(newdim[0],newdim[1])
+ tmp2(no_z_tmp)=1
+ tmp = tmp and tmp2
+ endfor
+ no_z = where(tmp eq 1)
+ Weight = total(Weight,3)
+ endif
+ ATM = dblarr(work_size)
+ ATM[no_z] = Weight[no_z]
+
+ ;;; initial value of x
+ xtmp = mean(flux)*dim[2]/total(Weight(no_z))
+ x = dblarr(work_size)
+ x[no_z] = xtmp
+ Wtmp = dblarr(work_size)
+ Wtmp(no_z) = dim[2]/Weight(no_z)
+ Weight = Wtmp
+endif else begin
+ work_index = transpose(lindgen(product(work_size)))
+ no_z = work_index
+endelse
+
+;;;;;;;;;;;;;;;;;;;;;
+; vector allocation ;
+;;;;;;;;;;;;;;;;;;;;;
+if errflag then begin
+ err = dblarr(maxit+1)
+ obj = obj/scaling
+ obj_sum = total(obj*obj)
+;help, obj, obj_sum
+endif
+Discr = dblarr(maxit+1)
+times = dblarr(maxit+1)
+
+;;;;;;;;;;;;;;;;
+; start of SGP ;
+;;;;;;;;;;;;;;;;
+; projection of the initial point
+x_nonpos = where(x lt 0, n_nonpos)
+if (n_nonpos gt 0) then x(x_nonpos) = 0.
+dummy=total(x,/double)
+x=x*pr_mask ;projection on user-defined mask
+x=x/total(x,/double)*dummy
+
+; error
+if errflag then begin
+ e = x(work_index) - obj[*]
+ err[0] = sqrt(total(e*e)/obj_sum)
+;help, e, err[0]
+endif
+
+if bound then begin
+ tmp = dblarr([work_size,dim[2]])
+ for i = 0, dim[2]-1 do begin
+ tmp[*,*,i] = Mask
+ endfor
+ multi_index = where(tmp gt 0)
+ tmp = make_array(work_size[0],work_size[1],dim[2])
+ tmp[multi_index] = gn[*]
+ gn = tmp
+ tmp = make_array(work_size[0],work_size[1],dim[2],/double,value=median(bg))
+ tmp[multi_index] = bg[*]
+ bg = tmp
+endif else begin
+ multi_index = lindgen(dim[2]*product(work_size))
+endelse
+
+; objective function value
+x_tf = afunction(x,TF)
+den = x_tf + bg
+temp = gn/den
+g = afunction(temp,CTF)
+
+if dim[2] gt 1 then begin
+ g = total(g,3)
+endif
+if bound then begin
+ g = ATM - g
+endif else begin
+ g = ONE - g
+endelse
+fv = total(gn(multi_index)*alog(temp(multi_index))) + total(x_tf(multi_index)) - total(flux)
+
+;;;;; bounds for the scaling matrices ;;;;;;;;
+temp = afunction(gn,CTF)
+tmp = temp[*,*,0]
+bg_tmp = bg[*,*,0]
+y = (flux[0]/(flux[0] + total(bg_tmp[work_index]) ))*tmp(work_index)
+y_min = min(y[where(y gt 0)])
+y_max = max(y)
+for i = 1, dim[2]-1 do begin
+ tmp = temp[*,*,i]
+ bg_tmp = bg[*,*,i]
+ y = (flux[i]/(flux[i] + total(bg_tmp[work_index])))*tmp(work_index)
+ y_min = min([min(y[where(y gt 0)]),y_min])
+ y_max = max([max(y),y_max])
+endfor
+
+X_low_bound = y_min ; Lower bound for the scaling matrix
+X_upp_bound = y_max ; Upper bound for the scaling matrix
+
+if X_upp_bound/X_low_bound lt 50. then begin
+ X_low_bound = X_low_bound/10.
+ X_upp_bound = X_upp_bound*10.
+endif
+; discrepancy
+Discr[0] = Discr_coeff * fv
+
+; scaling matrix
+if initflag eq 0 then begin
+ XX = make_array(size(x,/dimensions), value = 1., /double)
+endif else begin
+ XX = x
+ ; bounds
+ XX_toolow = where(XX lt X_low_bound, nlow)
+ if (nlow gt 0) then XX(XX_toolow) = X_low_bound
+ XX_toohigh = where(XX gt X_upp_bound, nhigh)
+ if (nhigh gt 0) then XX(XX_toohigh) = X_upp_bound
+ if bound then XX = Weight*XX
+end
+
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+; tolerance for stop criterion
+;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+
+loop = 1B
+while loop do begin
+
+ Valpha[0:Malpha-2] = Valpha[1:Malpha-1]
+ if M gt 1 then Fold[0:M-2] = Fold[1:M-1]
+ Fold[M-1] = fv
+ ; Step 2.1
+ y = x - alpha*XX*g
+ ; projection
+ y_nonpos = where(y lt 0, n_nonpos)
+ if (n_nonpos gt 0) then y(y_nonpos) = 0
+ y=y*pr_mask ;projection on user-defined mask
+ d = y - x
+ gd = total(d*g)
+ lam = 1.
+ ; Step 2.2
+ fcontinue = 1B
+ ; exploiting linearity
+ d_tf = afunction(d,TF)
+ fr = max(Fold)
+
+ while fcontinue do begin
+ xplus = x + lam*d
+ x_tf_try = x_tf + lam*d_tf
+ den = x_tf_try + bg
+ temp = gn/den
+ fv = total( gn(multi_index)* alog(temp(multi_index)) ) + total( x_tf_try(multi_index) ) - total(flux)
+ ; Step 2.3
+ if ( fv le fr + beta * lam * gd OR lam lt 1.e-12) then begin
+ x = xplus ;clear xplus
+ sk = lam*d
+ x_tf = x_tf_try
+ gtemp = afunction(temp,CTF)
+ if dim[2] gt 1 then begin
+ gtemp = total(gtemp,3)
+ endif
+ if bound then begin
+ gtemp = ATM - gtemp
+ endif else begin
+ gtemp = ONE - gtemp
+ endelse
+
+ yk = gtemp - g
+ g = gtemp
+ fcontinue = 0B
+ endif else begin
+ lam = lam * theta
+ endelse
+ endwhile
+
+ if (fv ge fr AND verb gt 0) then print, 'WARNING: fv >= fr'
+
+ ; Step 3
+ XX = x
+ ; bounds
+ XX_toolow = where(XX lt X_low_bound, nlow)
+ if (nlow gt 0) then XX(XX_toolow) = X_low_bound
+ XX_toohigh = where(XX gt X_upp_bound, nhigh)
+ if (nhigh gt 0) then XX(XX_toohigh) = X_upp_bound
+ if bound then XX = Weight*XX
+ D = dblarr(work_size)
+ D[no_z] = 1./XX[no_z]
+
+ sk2 = sk*D
+ yk2 = yk*XX
+
+ bk = total(sk2*yk)
+ ck = total(yk2*sk)
+
+ if (bk le 0) then begin
+ alpha1 = min([10.*alpha,alpha_max])
+ endif else begin
+ alpha1BB = total(sk2*sk2)/bk
+ alpha1 = min([alpha_max, max([alpha_min, alpha1BB])])
+ end
+ if (ck le 0) then begin
+ alpha2 = min([10.*alpha,alpha_max])
+ endif else begin
+ alpha2BB = ck/total(yk2*yk2)
+ alpha2 = min([alpha_max, max([alpha_min, alpha2BB])])
+ end
+
+ Valpha[Malpha-1] = alpha2
+
+ if (iter le 20) then begin
+ alpha = min(Valpha)
+ endif else begin
+ if (alpha2/alpha1 lt tau) then begin
+ alpha = min(Valpha)
+ tau = tau*0.9
+ endif else begin
+ alpha = alpha1
+ tau = tau*1.1
+ endelse
+ endelse
+ times[iter] = systime(2) - t0
+ iter = iter + 1
+
+ alpha = double(float(alpha))
+
+ if errflag then begin
+ e = x(work_index) - obj[*]
+ err[iter-1] = sqrt(total(e*e)/obj_sum)
+;help, e, err[iter-1]
+ endif
+
+ Discr[iter-1] = Discr_coeff * fv
+;;;;;;;;;;;;;;;;;
+; stop criteria ;
+;;;;;;;;;;;;;;;;;
+
+ case stopcrit of
+ 1: begin
+ if verb gt 0 then begin
+ print, 'it ', iter-1, ' of ', maxit
+ endif
+ end
+ 2: begin
+ normstep = total(sk*sk)
+ loop = (normstep gt tol*total(x*x))
+ if verb gt 0 then begin
+ print, 'it ', iter-1, ', || x_k - x_k-1 || ^2 / || x_k || ^2 ', normstep, ', tol ', tol
+ endif
+ end
+ 3: begin
+ reldecrease = abs(fv - Fold[M-1])/abs(fv)
+ loop = (reldecrease gt tol)
+ if verb gt 0 then begin
+ print, 'it ', iter-1, ', | f_k - f_k-1 | / | f_k | ', reldecrease, ', tol ', tol
+ endif
+ end
+ 4: begin
+ if iter ne 1 then begin
+ loop = (Discr[iter-1] gt tol)
+ endif else begin
+ loop = 1
+ endelse
+ if verb gt 0 then begin
+ print, 'it ', iter-1, ', D_k ', Discr[iter-1], ', tol ', tol
+ endif
+ end
+ endcase
+
+ if iter gt maxit then loop = 0B
+
+ if verb gt 0 then begin
+ print, 'Iteration:', iter, ' Fobj:', fv, ' Alpha:', alpha, ' Lambda:', lam, ' Discr:', Discr[iter-1]
+ endif
+
+endwhile
+;code by ALC - 2013/08/09:
+;compute residual and its statistics
+tmp = afunction(x,TF)+ bg
+res = (gn - tmp) / sqrt(tmp)*sqrt(scaling)
+res = res(work_index)
+res = reform(res,obj_size)
+histogauss, res, res_stat, /NOPLOT
+
+;reconstructed object:
+x = x(work_index)
+x = reform(x,obj_size)
+x = x * scaling
+
+if errflag then err = err[0:iter-1]
+
+Discr = Discr[0:iter-1]
+times = times[0:iter-1]
+iter = iter - 1
+
+END
+
+; ==============================================================================
+; End of SGP_deblurring.pro file - IRMA package
+; ==============================================================================
diff --git a/Patch_lib/utility_for_noi_table.pro b/Patch_lib/utility_for_noi_table.pro
new file mode 100644
index 0000000000000000000000000000000000000000..bc9fae236e65edc15ff8ddfce2c8047e0ccda369
--- /dev/null
+++ b/Patch_lib/utility_for_noi_table.pro
@@ -0,0 +1,76 @@
+;+
+;RESIZE A TABLE OF NOI (NUMBER OF ITERATION)
+;WITH INTERPOLATION FROM N to N'
+;(N > N')
+;
+; USAGE:
+; output = utility_for_noi_table(Nx, Ny, br,bl,ul,ur, Mx, My)
+;
+; INPUT:
+; Nx = number of cols of input table
+; Ny = number of rows " "
+; br = bottom right value " "
+; bl = bottom left value " "
+; ur = upper right value " "
+; ul = upper left value " "
+; Mx = number of cols of output table
+; My = number of rows of output table
+;
+; OUTPUT:
+; a Mx X My array containing the computed values
+;
+; KEYWORDS:
+; VERBOSE : print every phase of the computation (input table and output table)
+;
+; EXAMPLE:
+; my_table = utility_for_noi_table(15,15,10,20,20,50,9,9,/VERBOSE)
+;
+;
+;-
+
+pro my_print, tab, Nx, Ny
+;print an 1D array as a Nx X Ny table
+;rotate(reform(tab,Nx,Ny),7)
+for j=0,Ny-1 do begin ;;ROWS
+ tmp_string = ""
+ for i=0,Nx-1 do begin
+ tmp_string=tmp_string+strtrim(tab[i+Nx*(Ny-j-1)],1)+" "
+ endfor
+ print, tmp_string
+endfor
+end
+
+
+
+function utility_for_noi_table, Nx, Ny, bl,br,ul,ur, Mx, My, VERBOSE=verbose
+
+if KEYWORD_SET(verbose) then verb = 1B else verb = 0B
+
+;;;;;;;;;;;;;;;;;
+;input section
+;;;;;;;;;;;;;;;;;
+
+;definition of the table
+N=Nx*Ny
+a=[[ul,ur],[bl,br]] ;as done within "tecno_inaf"
+a=float(rotate(a,7)) ; as done within "tecno_inaf"
+tab=make_array(N,/LONG)
+for k=0,N-1 do begin
+ i=(k mod Nx)
+ j=(k / Nx)
+ d1=a[0,0]+(a[1,0]-a[0,0])/(Nx-1)*i
+ d2=a[0,1]+(a[1,1]-a[0,1])/(Nx-1)*i
+ d3=round(d1+(d2-d1)/(Ny-1)*j)
+ tab[k]=d3
+endfor
+
+;print, tab[*]
+if (verb) then my_print, tab, Nx, Ny
+if (verb) then print, " "
+
+tab2=rotate(congrid(rotate(reform(tab,Nx,Ny),7),Mx,My,/INTERP),7)
+;print, tab2[*]
+if (verb) then my_print, tab2, Mx, My
+
+return, tab2
+end
diff --git a/Patch_lib/website.html b/Patch_lib/website.html
new file mode 100644
index 0000000000000000000000000000000000000000..002a542b9b0de5c764687232c4a24a49414f23c5
--- /dev/null
+++ b/Patch_lib/website.html
@@ -0,0 +1,9 @@
+<html>
+<head>
+<meta http-equiv="Content-Type" content="text/html; charset=ISO-8859-1">
+<title>TECNO INAF website</title>
+<meta HTTP-EQUIV="refresh" content="0; url=http://www.airyproject.eu/Tutorials/Patch/">
+</head>
+<body>
+</body>
+</html>
diff --git a/Patch_lib/zero_padding.pro b/Patch_lib/zero_padding.pro
new file mode 100644
index 0000000000000000000000000000000000000000..785c161ac4f01b0fa70ee79f76042784f078dfb5
--- /dev/null
+++ b/Patch_lib/zero_padding.pro
@@ -0,0 +1,68 @@
+; $Id: zero_padding last version: 2012/01/20 Andrea La Camera $
+;
+;+
+; NAME:
+; zero_padding
+;
+; PURPOSE:
+; To enlarge the size of a 2d/3d array by using zero_padding. The final dimensions
+; are dim1 x dim2. 3d arrays are padded slice per slice, i.e. the third dimension
+; will not be changed.
+;
+; CATEGORY:
+; utility
+;
+; CALLING SEQUENCE:
+; result = zero_padding(img, dim1, dim2, [value=value])
+;
+; INPUT:
+; img : initial image, can be either 2D or 3D (cube of images)
+; dim1,dim2 : output dimensions (dim1xdim2)
+;
+; OPTIONAL INPUT:
+; value : value of the extended boundary (default = 0)
+;
+; OUTPUT:
+; result : the zero-padded image(s)
+;
+; ROUTINE MODIFICATION HISTORY:
+; routine written: May 2007,
+; Andrea La Camera (DISI) [lacamera@disi.unige.it]
+; modifications : November 2011,
+; Andrea La Camera (DISI) [lacamera@disi.unige.it]
+; -a few modifications
+; January 2012,
+; Andrea La Camera (DISI) [lacamera@disi.unige.it]
+; -added option "value"
+;-
+
+
+function zero_padding, img, dim1, dim2, value=value
+
+IF N_PARAMS() NE 3 then begin
+ print, "Syntax:"
+ print, " result = zero_padding(image, N, M)"
+ print, "where image is the input; N and M are the dimensions of the output"
+ return, -1
+endif
+on_error,2
+
+IF KEYWORD_SET(value) THEN ext_value=value ELSE ext_value=0
+
+N=(size(img))[1]
+M=(size(img))[2]
+if ((size(img))[0] EQ 1) OR ((size(img))[0] GT 3) then message, "Only 2d/3d arrays"
+if (dim1 LT N) OR (dim2 LT M) then message, $
+ "Final dimensions must be greater than or equal to the initial dimensions."
+if (size(img))[0] EQ 2 then begin
+ img_zp=make_array(dim1,dim2,type=(size(img,/type)), value=ext_value)
+ img_zp[floor((dim1-N)/2.),floor((dim2-M)/2.)]=img
+endif else begin
+ P=(size(img))[3]
+ img_zp=make_array(dim1,dim2,P,type=(size(img,/type)), value=ext_value)
+ for j=0,P-1 do img_zp[floor((dim1-N)/2.),floor((dim2-M)/2.),j]=img[*,*,j]
+endelse
+
+
+return, img_zp
+end
diff --git a/Patch_setup/patch_env.sh b/Patch_setup/patch_env.sh
new file mode 100644
index 0000000000000000000000000000000000000000..c33fcb705e87e4be96b5d81d65b927cd0ab90fda
--- /dev/null
+++ b/Patch_setup/patch_env.sh
@@ -0,0 +1,3 @@
+#REPLACE THE STRING WITH YOUR INSTALLATION PATH
+# startup file path definition
+export IDL_STARTUP="/Users/andrea/Dropbox/Patch/Patch_setup/patch_startup.pro"
diff --git a/Patch_setup/patch_startup.pro b/Patch_setup/patch_startup.pro
new file mode 100644
index 0000000000000000000000000000000000000000..7bf7824bb211037bbcea2c60b4b973b0cd14545a
--- /dev/null
+++ b/Patch_setup/patch_startup.pro
@@ -0,0 +1,18 @@
+;;; PATCH path expansion
+!PATH=!PATH+':'+EXPAND_PATH('+/Users/andrea/Dropbox/Patch/')
+!HELP_PATH=!HELP_PATH+':'+EXPAND_PATH('+/Users/andrea/Dropbox/Patch/')
+
+;;; Required lib
+;;; 1) "The IDL Astronomy User's Library" downloadable from
+;;; http://idlastro.gsfc.nasa.gov
+!PATH=!PATH+':'+EXPAND_PATH('+/Users/andrea/Dropbox/DRS_NEW/astron/')
+
+;;; Other optional libraries / user-defined working dir
+;;; e.g. !PATH=!PATH+':'+EXPAND_PATH('+/Users/andrea/myworkingdir/')
+
+;; some optional IDL options
+device, decomposed=0
+device, retain=2
+
+;;; prompt modification
+!prompt="PATCH > "
diff --git a/README.txt b/README.txt
new file mode 100644
index 0000000000000000000000000000000000000000..643c68d7d093b7a48d4b28a708ecdf04fc6f7c79
--- /dev/null
+++ b/README.txt
@@ -0,0 +1,91 @@
+ SOFTWARE PATCH
+ v. 0.9
+ written by Andrea La Camera
+ (DIBRIS, University of Genova, Italy)
+ October 2014
+
+
+
+Software PATCH is an IDL software developed for the deconvolution of
+mono-pupil images with space-variant point-spread functions (PSFs). The input
+image is decomposed in partially overlapping sub-domains, the size of which will
+ depend on the properties of the PSF.
+
+The reconstruction is obtained by applying a deconvolution method with boundary
+effects compensation to each sub-domain. The final object is formed by a mosaic
+of the partial reconstructed images, after the removal of the common parts of
+adjacent sub-domains.
+
+Two algorithms are available: the well-known Richardson-Lucy and
+the Scaled Gradient Projection (SGP) developed by S. Bonettini et al.
+
+More information and a tutorial can be found on
+http://www.airyproject.eu/Tutorials/Patch/
+
+
+
+
+-----------------------------
+REQUIREMENTS OF THE SOFTWARE:
+-----------------------------
+- IDL version 7.0 (tested on 7.0, not guaranteed for IDL 8.0)
+- IDL Astronomy User's Library (downloadable from
+ http://idlastro.gsfc.nasa.gov/homepage.html)
+
+--------------------------
+INSTALLATION INSTRUCTIONS:
+--------------------------
+
+1 - Extract the archive "patch.zip" somewhere in your machine; for
+ instance, you can extract on your desktop. Depending on your OS, you will have:
+ - C:\Documents and Settings\andrea\Desktop\
+ - /Users/andrea/Desktop/
+ - /home/andrea/Desktop/
+
+ The main directory "PATCH" contains the main program
+ ("patch.pro"), a sub-folder for your data and two sub-folders
+ containing the library and the setup files. The main folder is organized as
+ follows:
+ [path_on_your_OS]/Patch
+ | - Patch_data
+ | - Patch_lib
+ | - Patch_setup
+ patch.pro
+ README.txt (this file)
+
+2a - UNIX/LINUX/MACOsX case:
+ i) Open the subfolder "Patch_setup" and edit
+ the "patch_env.sh", fixing the path of the idl_startup file.
+ Then edit the "patch_startup.pro" file and fix the path of your
+ working directories (PATCH, astrolib and, if you like, other folders
+ such as your Desktop, a tmp folder, etc.).
+ ii) Open your terminal and type:
+ $ source your_path_where_patch_is/Patch_setup/patch_env.sh
+iii) launch IDL. Note that every time you want to launch the "PATCH"
+ program from a new terminal session, you have to perform both step ii) and
+ iii).
+
+[Suggested] Edit your .bashrc (or your shell's configuration file) and add an
+ alias for PATCH. For example, you can create an alias called "patch":
+ alias patch='source /home/andrea/patch/patch_setup/patch_env.sh ; idl'
+ In such a way, you can skip both ii) and iii) steps and launch directly
+ Software PATCH by typing the "patch" command at your terminal prompt.
+
+
+2b - WINDOWS case (IDL Workbench case):
+ i) You can delete the "patch_env.sh" file.
+ ii) Then set the paths and the startup file in the IDL preferences
+ (Window->Preferences->IDL):
+ - add the "Patch" main folder to path and include subfolders;
+ - add the "astrolib" main folder to path and include subfolders;
+ - add the "Patch" main folder to help-path and include subfolders;
+ - set the "patch_startup.pro" file as IDL startup file.
+iii) Click on "Apply" and reset your IDL session (or exit and re-launch IDL).
+
+3 - Type "patch" at the IDL prompt to launch the Graphical User Interface
+ (GUI) of the program.
+
+
+---
+Last update: 2014-10-16
+Andrea La Camera (DIBRIS) [andrea.lacamera@unige.it]
diff --git a/patch.pro b/patch.pro
new file mode 100644
index 0000000000000000000000000000000000000000..0a92b091ed2b3a7e16e0ca273e988269bc0f55fb
--- /dev/null
+++ b/patch.pro
@@ -0,0 +1,93 @@
+;---patch.pro, v. 0.9, last modified: 2014-10-16, Andrea La Camera
+;+
+; NAME:
+; PATCH
+; PURPOSE:
+; PATCH program is an IDL software developed for the deconvolution of
+; mono-pupil images with space-variant point-spread functions (PSFs). The
+; input image is decomposed in partially overlapping sub-domains, the size
+; of which will depend on the properties of the PSF.
+; The reconstruction is obtained by applying a deconvolution method with
+; boundary effects compensation to each sub-domain. The final object is
+; formed by a mosaic of the partial reconstructed images, after the
+; removal of the common parts of adjacent sub-domains.
+; Two algorithms are available: the well-known Richardson-Lucy and
+; the Scaled Gradient Projection (SGP) developed by S. Bonettini et al.
+; More information and a tutorial can be found on
+; http://www.airyproject.eu/Tutorials/Patch/
+;
+; CALLING SEQUENCE:
+; 1) In order to start the Graphical User Interface (GUI) of the program
+; and to set all parameters, please enter at IDL prompt:
+; PATCH > patch
+; 2) In order to directly launch the deconvolution step, enter at IDL prompt:
+; PATCH > patch, path_param_file
+; where 'path_param_file' is the path of the param_file.sav that
+; contains all parameters necessary for the reconstruction (previously
+; saved from the GUI).
+;
+; INPUTS:
+; mode 1) --> no input
+; mode 2) --> path_param_file in the form of "/home/andrea/Desktop/mypar.sav"
+;
+; OPTIONAL INPUT KEYWORD:
+; None
+;
+; OUTPUTS:
+; None
+;
+; RESTRICTIONS:
+; None
+;
+; MODIFICATION HISTORY:
+; WRITTEN, Andrea La Camera (DIBRIS, Univ. of Genova), November 2011
+; v0.81 - Minor modifications in the GUI (2014-09-29, ALC)
+; v0.9 - Overlap management changed. Added a new routine for
+; this goal. Changes in patch_gui.pro and patch_deconv.pro (2014-10-16 ALC)
+;-
+
+pro patch, param_file
+
+on_error, 2
+COMMON image, img, header, N, M, a_dim, Nori, Mori
+COMMON background, bg ;bg must be N x M, as the input image "img" is.
+COMMON output, out, res_array, residual ; reconstructed object with original dimension Nori x Mori
+
+if keyword_set(param_file) then begin
+ ;message, "NOT YET IMPLEMENTED"
+ ;load the parameters (param)
+ restore, param_file, /VERBOSE
+ ;define the structure
+ state={par:param}
+ ; load image and define img, N, M
+ a=patch_load_image(state, 0, VERBOSE=0)
+ if a EQ 0 then message, 'ERROR on loading the input image file.'
+ ;check if zero_padding is neccessary...
+ if (state.par.cols*state.par.domx GT N) OR $
+ (state.par.rows*state.par.domy GT M) then begin
+ N=state.par.cols*state.par.domx
+ M=state.par.rows*state.par.domy
+ img = zero_padding(img, N, M)
+ endif
+ ;load (or define) the sky background array
+ if state.par.bg_choice EQ 0 then begin
+ bg = make_array(N,M,value=state.par.bg)
+ endif else begin
+ a=patch_load_background(state,0,VERBOSE=0)
+ if a EQ 0 then message, 'ERROR on loading the background file.'
+ endelse
+
+ ;RUN DECONVOLUTION!
+ if n_elements(img) EQ 0 then message, 'Input image undefined!'
+ if n_elements(bg) EQ 0 then message, 'Background undefined!'
+ t0=systime(1)
+ print, "Please, wait..."
+ patch_deconv, state.par, /SAVE_OUT
+ print, "...done!"
+ print, "Total deconvolution time (s) = ", systime(1)-t0
+ print, "Average time for each sub-domain = ", $
+ (systime(1)-t0)/(state.par.rows*state.par.cols)
+
+endif else patch_gui
+
+end