Commit 70a469db authored by Laura Schreiber's avatar Laura Schreiber
Browse files

First commit

parent a920aaff
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
; ==============================================================================
; 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
;+
; 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 ; ----------------------------------------------------------------------------
;+
; NAME:
; FPUFIX
;
; PURPOSE:
;
; This is a utility routine to examine a variable and fix problems
; that will create floating point underflow errors.