pax_global_header 0000666 0000000 0000000 00000000064 14220535751 0014516 g ustar 00root root 0000000 0000000 52 comment=caa897d8b7e49912d385035411f3169c87082379 patch-caa897d8b7e49912d385035411f3169c87082379/ 0000775 0000000 0000000 00000000000 14220535751 0017210 5 ustar 00root root 0000000 0000000 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/ 0000775 0000000 0000000 00000000000 14220535751 0021075 5 ustar 00root root 0000000 0000000 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/afunction.pro 0000664 0000000 0000000 00000007331 14220535751 0023611 0 ustar 00root root 0000000 0000000 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 ; ============================================================================== patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/coyote_library/ 0000775 0000000 0000000 00000000000 14220535751 0024123 5 ustar 00root root 0000000 0000000 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/coyote_library/cgscalevector.pro 0000664 0000000 0000000 00000021127 14220535751 0027474 0 ustar 00root root 0000000 0000000 ; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/coyote_library/error_message.pro 0000664 0000000 0000000 00000030506 14220535751 0027506 0 ustar 00root root 0000000 0000000 ;+ ; 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 ; ---------------------------------------------------------------------------- patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/coyote_library/fpufix.pro 0000664 0000000 0000000 00000012014 14220535751 0026144 0 ustar 00root root 0000000 0000000 ;+ ; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/coyote_library/gmascl.pro 0000664 0000000 0000000 00000016454 14220535751 0026125 0 ustar 00root root 0000000 0000000 ;+ ; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/coyote_library/logscl.pro 0000664 0000000 0000000 00000016272 14220535751 0026140 0 ustar 00root root 0000000 0000000 ;+ ; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/coyote_library/scale_vector.pro 0000664 0000000 0000000 00000021242 14220535751 0027317 0 ustar 00root root 0000000 0000000 ; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_check_overlap.pro 0000664 0000000 0000000 00000004006 14220535751 0025603 0 ustar 00root root 0000000 0000000 ;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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_deconv.pro 0000664 0000000 0000000 00000025025 14220535751 0024260 0 ustar 00root root 0000000 0000000 ;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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_deconv_rl.pro 0000664 0000000 0000000 00000001055 14220535751 0024752 0 ustar 00root root 0000000 0000000 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_deconv_sgp.pro 0000664 0000000 0000000 00000000742 14220535751 0025130 0 ustar 00root root 0000000 0000000 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_gui.pro 0000664 0000000 0000000 00000105124 14220535751 0023565 0 ustar 00root root 0000000 0000000 ;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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_img_define.pro 0000664 0000000 0000000 00000013317 14220535751 0025071 0 ustar 00root root 0000000 0000000 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_img_options.pro 0000664 0000000 0000000 00000005531 14220535751 0025331 0 ustar 00root root 0000000 0000000 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_load_background.pro 0000664 0000000 0000000 00000004122 14220535751 0026113 0 ustar 00root root 0000000 0000000 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_load_image.pro 0000664 0000000 0000000 00000002441 14220535751 0025060 0 ustar 00root root 0000000 0000000 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_load_psfs.pro 0000664 0000000 0000000 00000007225 14220535751 0024756 0 ustar 00root root 0000000 0000000 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_noi_table.pro 0000664 0000000 0000000 00000005203 14220535751 0024732 0 ustar 00root root 0000000 0000000 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_noi_table_compute.pro 0000664 0000000 0000000 00000001640 14220535751 0026467 0 ustar 00root root 0000000 0000000 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_psf_ee.pro 0000664 0000000 0000000 00000001654 14220535751 0024245 0 ustar 00root root 0000000 0000000 ;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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_set_gui.pro 0000664 0000000 0000000 00000005270 14220535751 0024441 0 ustar 00root root 0000000 0000000 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_show_image.pro 0000664 0000000 0000000 00000006443 14220535751 0025127 0 ustar 00root root 0000000 0000000 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_show_out.pro 0000664 0000000 0000000 00000003703 14220535751 0024650 0 ustar 00root root 0000000 0000000 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/patch_show_res.pro 0000664 0000000 0000000 00000003731 14220535751 0024633 0 ustar 00root root 0000000 0000000 ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; 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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/sgp_deblurring.pro 0000664 0000000 0000000 00000046621 14220535751 0024636 0 ustar 00root root 0000000 0000000 ;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 ; ============================================================================== patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/utility_for_noi_table.pro 0000664 0000000 0000000 00000003316 14220535751 0026207 0 ustar 00root root 0000000 0000000 ;+ ;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 patch-caa897d8b7e49912d385035411f3169c87082379/Patch_lib/website.html 0000664 0000000 0000000 00000000360 14220535751 0023424 0 ustar 00root root 0000000 0000000