;+ ; NAME: ; Sharpen ; ; PURPOSE: ; ; This function sharpens an image using a Laplacian kernel. ; The final result is color adjusted to match the histogram ; of the input image. ; ; 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: ; ; Image Processing ; ; CALLING SEQUENCE: ; ; sharp_image = Sharpen(image) ; ; INPUTS: ; ; image - The input image to be sharpened. Assumed to be a 2D byte array. ; ; OUTPUTS: ; ; sharp_image - The sharpened image. ; ; INPUT KEYWORDS: ; ; KERNEL -- By default the image is convolved with this 3-by-3 Laplacian kernel: ; [ [-1, -1, -1], [-1, +8, -1], [-1, -1, -1] ]. You can pass in any kernel ; of odd width. The filtered image is added back to the original image to provide ; the sharpening effect. ; ; DISPLAY -- If this keyword is set a window is opened and the details of the sharpening ; process are displayed. ; ; OUTPUT KEYWORDS: ; ; None. ; ; DEPENDENCIES: ; ; None. ; ; METHOD: ; ; This function is based on the Laplacian kernel sharpening method on pages 128-131 ; of Digital Image Processing, 2nd Edition, Rafael C. Gonzalez and Richard E. Woods, ; ISBN 0-20-118075-8. ; ; EXAMPLE: ; ; There is an example program at the end of this file. ; ; MODIFICATION HISTORY: ; ; Written by David W. Fanning, January 2003. ; Updated slightly to use Coyote Library routines. 3 Dec. 2010. DWF. ; Modified the example to work with cgImage. 29 March 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 Sharpen_HistoMatch, image, histogram_to_match ; Error handling. Catch, theError IF theError NE 0 THEN BEGIN Catch, /Cancel ; Get the call stack and the calling routine's name. Help, Calls=callStack IF Float(!Version.Release) GE 5.2 THEN $ callingRoutine = (StrSplit(StrCompress(callStack[1])," ", /Extract))[0] ELSE $ callingRoutine = (Str_Sep(StrCompress(callStack[1])," "))[0] ; Print a traceback. Help, /Last_Message, Output=traceback Print,'' Print, 'Traceback Report from ' + StrUpCase(callingRoutine) + ':' Print, '' FOR j=0,N_Elements(traceback)-1 DO Print, " " + traceback[j] IF N_Elements(image) NE 0 THEN RETURN, image ELSE RETURN, -1L ENDIF ; We require two input parameters. IF N_Params() NE 2 THEN Message, 'Two arguments required. Please read the program documentation.' ; Must have 2D image array. IF Size(image, /N_Dimensions) NE 2 THEN Message, 'Image argument must be 2D. Returning.' ; Is the histogram_to_match variable a 1D or 2D array? Branch accordingly. CASE Size(histogram_to_match, /N_Dimensions) OF 1: BEGIN IF N_Elements(histogram_to_match) NE 256 THEN $ Message, 'Histogram to match has incorrect size. Returning.' match_histogram = histogram_to_match END 2: match_histogram = Histogram(Byte(histogram_to_match), Min=0, Max=255, Binsize=1) ELSE: Message, 'Histogram to match has incorrect number of dimensions. Returning.' ENDCASE ; Calculate the histogram of the input image. h = Histogram(Byte(image), Binsize=1, Min=0, Max=255) ; Make sure the two histograms have the same number of pixels. This will ; be a problem if the two images are different sizes, you are matching a ; histogram from an image subset, etc. totalPixels = Float(N_Elements(image)) totalHistogramPixels = Float(Total(match_histogram)) IF totalPixels NE totalHistogramPixels THEN $ factor = totalPixels / totalHistogramPixels ELSE $ factor = 1.0 match_histogram = match_histogram * factor ; Find a mapping from the input pixels to the transformation function s. s = FltArr(256) FOR k=0,255 DO BEGIN s[k] = Total(h(0:k) / totalPixels) ENDFOR ; Find a mapping from input histogram to the transformation function v. v = FltArr(256) FOR q=0,255 DO BEGIN v[q] = Total(match_histogram(0:q) / Total(match_histogram)) ENDFOR ; Find probablitly density function z from v and s. z = BytArr(256) FOR j=0,255 DO BEGIN i = Where(v LT s[j], count) IF count GT 0 THEN z[j] = (Reverse(i))[0] ELSE z[j]=0 ENDFOR ; Create the matched image. matchedImage = z[Byte(image)] RETURN, matchedImage END ; ---------------------------------------------------------------------------- FUNCTION Sharpen, image, Display=display, Kernel=kernel ; Error handling. Catch, theError IF theError NE 0 THEN BEGIN Catch, /Cancel ; Get the call stack and the calling routine's name. Help, Calls=callStack IF Float(!Version.Release) GE 5.2 THEN $ callingRoutine = (StrSplit(StrCompress(callStack[1])," ", /Extract))[0] ELSE $ callingRoutine = (Str_Sep(StrCompress(callStack[1])," "))[0] ; Print a traceback. Help, /Last_Message, Output=traceback Print,'' Print, 'Traceback Report from ' + StrUpCase(callingRoutine) + ':' Print, '' FOR j=0,N_Elements(traceback)-1 DO Print, " " + traceback[j] IF N_Elements(image) NE 0 THEN RETURN, image ELSE RETURN, -1L ENDIF ; If an image is not provided. Issue an error message. IF N_Elements(image) EQ 0 THEN $ Message, 'A 2D image is required as an argument.' IF Size(image, /N_Dimensions) NE 2 THEN Message, 'Image must be a 2D array in this program.' ; Resize the image, if required. previewSize = 512 wxsize = previewSize wysize = previewSize ; Set up the convolution kernel for Laplacian filtering. IF N_Elements(kernel) EQ 0 THEN BEGIN k = Replicate(-1, 3, 3) k[1,1] = 8 ENDIF ELSE BEGIN s = Size(kernel, /Dimensions) IF s[0] MOD 2 NE 1 THEN Message, 'Kernel must be an odd width.' k = kernel ENDELSE ; Are we doing a display? IF Keyword_Set(display) THEN BEGIN s = Size(image, /Dimensions) xsize = s[0] ysize = s[1] needresize = 1 IF xsize NE ysize THEN BEGIN needresize = 1 aspect = Float(ysize) / xsize IF aspect LT 1 THEN BEGIN wxsize = previewSize wysize = (previewSize * aspect) < previewSize ENDIF ELSE BEGIN wysize = previewSize wxsize = (previewSize / aspect) < previewSize ENDELSE ENDIF Window, /Free, XSize=2*wxsize, YSize=2*wysize, Title='Image Sharpening-Laplacian' ENDIF ELSE needresize = 0 ; Need a resize? IF needresize THEN thisImage = Byte(Congrid(image, wxsize, wysize)) ELSE $ thisImage = image ; Display the original image. IF Keyword_Set(display) THEN BEGIN $ cgImage, thisImage, 0, 0, /TV XYOUTS, wxsize/2, 10, /Device, 'Original Image', Font=0, $ Alignment=0.5, Color=cgColor('red6') ENDIF ; Create the Laplacian filtered image. filteredImage = Convol(Float(thisImage), k, Center=1, /Edge_Truncate, /NAN) ; Display the filtered image. IF Keyword_Set(display) THEN BEGIN fimage = Convol(thisImage, k, Center=1, /Edge_Truncate, /NAN) cgImage, fimage, wxsize, wysize, /TV XYOUTS, (2*wxsize/4)*3, wysize + 10, /Device, 'Filtered Image', Font=0, $ Alignment=0.5, Color=cgColor('red6') ENDIF ; Scale the Laplacian filtered image. Note conversion of ; image to integer values and use of 255 as a FLOAT value. filteredImage = filteredImage - (Min(filteredImage)) filteredImage = filteredImage * (255./Max(filteredImage)) IF Keyword_Set(display) THEN BEGIN cgImage, filteredImage, 0, wysize, /TV XYOUTS, wxsize/2, wysize + 10, /Device, 'Scaled Filter', Font=0, $ Alignment=0.5, Color=cgColor('red6') ENDIF ; Create the sharpened image by adding the Laplacian filtered image ; back to the original image and re-scaling. sharpened = thisImage + filteredImage sharpened = sharpened - (Min(sharpened)) sharpened = sharpened * (255./Max(sharpened)) ; Adjust the sharpened image to match the histogram of the original. adjusted = Sharpen_HistoMatch(sharpened, image) ; Display the adjusted image. IF Keyword_Set(display) THEN BEGIN cgImage, BytScl(adjusted), wxsize, 0, /TV XYOUTS, (2*wxsize/4)*3, 10, /Device, 'Sharpened Image', Font=0, $ Alignment=0.5, Color=cgColor('red6') ENDIF RETURN, adjusted END PRO Example image = cgDemoData(13) s = Size(image, /Dimensions) LoadCT, 0, /Silent Window, /Free, XSize=s[0]*2, YSize=s[1], Title='Image Sharpening' cgImage, image, 0, /TV XYOuts, 0.25, 0.1, /Normal, Alignment=0.5, 'Original Image', Font=0, Color=cgColor('red6') cgImage, Sharpen(image), 1, /NoErase, /TV XYOuts, 0.75, 0.1, /Normal, Alignment=0.5, 'Sharpened Image', Font=0, Color=cgColor('red6') END