home *** CD-ROM | disk | FTP | other *** search
Text File | 1997-07-08 | 43.6 KB | 1,364 lines |
- ;$Id: d_reconstr.pro,v 1.25 1997/04/18 23:01:28 tremblay Exp $
- ;
- ; Copyright (c) 1997, Research Systems, Inc. All rights reserved.
- ; Unauthorized reproduction prohibited.
- ;
- ;+
- ; FILE:
- ; d_reconstr.pro
- ;
- ; CALLING SEQUENCE: d_reconstr
- ;
- ; PURPOSE:
- ; Shows reconstruction techniques for images.
- ; (Computerized tomography)
- ;
- ; MAJOR TOPICS: Visualization and data analysis
- ;
- ; CATEGORY:
- ; IDL 5.0
- ;
- ; INTERNAL FUNCTIONS and PROCEDURES:
- ; pro DemoMenuChoice - Identify the menu button choice
- ; pro DemoMenuCreate - Create the menu bar
- ; pro None - Set the filter to none
- ; pro Ramlak - Compute the Ramlak filter
- ; pro Shepp_Logan - Compute the Shepp_Logan filter
- ; pro Lp_Cosine - Compute the Lp_Cosine filter
- ; pro Gen_Hamming - Compute the Gen_Hamming filter
- ; pro DEllipse - Draw an ellipse
- ; pro DCircle - Draw a circle
- ; pro DPoly - Draw a polygon
- ; pro Make_Phantom - Create a phantom object
- ; pro Reconstruct_It - Do the image reconstruction
- ; pro Recon_Redraw - Redraw the windows
- ; pro Rec_Demo_Event - Event handler
- ; pro Reconstr_Cleanup - Cleanup
- ; pro D_Reconstr - Main procedure
- ;
- ; EXTERNAL FUNCTIONS, PROCEDURES, and FILES:
- ; pro gettips - Read the tip file
- ; pro widtips - Create the tip widgets
- ; pro sizetips - Size the tip widgets
- ; pro sizetips - Change the tips text
- ; reconstr.txt
- ; reconstr.tip
- ; ctscan.dat
- ;
- ; REFERENCE: IDL Reference Guide, IDL User's Guide
- ;
- ; NAMED STRUCTURES:
- ; none.
- ;
- ; COMMON BLOCS:
- ; tomodemo_com
- ;
- ; MODIFICATION HISTORY:
- ; 1/94, DS - Written.
- ; 1/97, DAT - New GUI, IDL Style Guide.
- ;-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Given a uservalue from a menu button created
- ; by MenuCreate, the function returns the index
- ; of the choice within the category. Set the
- ; selected menu button to insensitive to signify
- ; selection, and set all other choices for the
- ; category to sensitive.
-
- function DemoMenuChoice, $
- Eventval, $ ; IN: uservalue from seleted menu button
- MenuItems, $ ; IN: menu item array, as returned by MenuCreate
- MenuButtons ; IN: button array as returned by MenuCreate
-
- ; Get the name less the last qualifier.
- ;
- i = STRPOS(eventval, '|', 0)
- while (i ge 0) do begin
- j = i
- i = STRPOS(eventval, '|', i+1)
- endwhile
-
- base = STRMID(eventval, 0, j+1) ;common buttons, includes last |
-
- ; Get the button sharing this basename.
- ;
- buttons = WHERE(STRPOS(MenuItems, base) eq 0)
-
- ; Get the index of selected item.
- ;
- this = (WHERE(eventval eq MenuItems))(0)
-
- ; For each button in category, sensitize.
- ;
- for i=0, N_ELEMENTS(buttons)-1 do begin
- index = buttons(i)
- WIDGET_CONTROL, MenuButtons(buttons(i)), SENSITIVE=index ne this
- endfor
-
- ; Return the Selected button's index.
- ;
- RETURN, this - buttons(0)
-
- end
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Create a menu from a string descriptor (MenuItems).
- ; Return the parsed menu items in MenuItems (overwritten),
- ; and the array of corresponding menu buttons in MenuButtons.
- ;
- ; MenuItems = (input/output), on input the menu structure in the form of
- ; a string array. Each button is an element, encoded as follows:
- ; Character 1 = integer bit flag. Bit 0 = 1 to denote a button with
- ; children. Bit 1 = 2 to denote this is the last child of its
- ; parent. Bit 2 = 4 to show that this button should initially
- ; be insensitive, to denote selection. Any combination of bits
- ; may be set.
- ; On RETURN, MenuItems contains the fully qualified button names.
- ; Characters 2-end = Menu button text. Text should NOT contain the character
- ; |, which is used to delimit menu names.
- ; MenuButtons = (output) button widget id's of the created menu.
- ; Bar_base = (input) ID of menu base.
- ; Prefix = prefix for this menu's button names. If omitted, no
- ; prefix.
- ;
- ;
- ; Example:
- ; MenuItems = ['1File', '0Save', '2Exit', $
- ; '1Edit', '3Cut', $
- ; '3Help']
- ; Creates a menu with three top level buttons (file, edit and help).
- ; File has 2 choices (save and exit), Edit has one choice, and help has none.
- ; On RETURN, MenuItems contains the fully qualified menu button names
- ; in a string array of the form: ['<Prefix>|File', '<Prefix>|File|Save',
- ; '<Prefix>|File|Exit', '<Prefix>|Edit',..., etc. ]
- ;
- pro DemoMenuCreate, $
- MenuItems, $ ; IN/OUT: menu structure/button names
- MenuButtons, $ ; OUT: button widget identifier
- Bar_base, $ ; IN: menu bar base identifier
- Prefix=prefix ; IN: (opt) prefix of the menu's button names
-
- ; Initialize working variables and arrays.
- ;
- level = 0
- parent = [bar_base, 0, 0, 0, 0, 0]
- names = STRARR(5)
- lflags = INTARR(5)
- MenuButtons = LONARR(N_ELEMENTS(MenuItems))
-
- if (N_ELEMENTS(prefix)) then begin
- names(0) = prefix + '|'
- endif else begin
- names(0) = '|'
- endelse
-
- for i = 0, N_ELEMENTS(MenuItems)-1 do begin
- flag = FIX(STRMID(MenuItems(i), 0, 1))
- txt = STRMID(MenuItems(i), 1, 100)
- uv = ''
-
- for j = 0, level do uv = uv + names(j)
-
- ; Set the fully qualified name in the menuItems array.
- ;
- MenuItems(i) = uv + txt
-
- ; Create the menu bar buttons.
- ;
- MenuButtons(i) = WIDGET_BUTTON(parent(level), $
- VALUE=txt, UVALUE=uv+txt, $
- MENU=flag and 1, HELP=txt eq 'About')
-
- if ((flag and 4) ne 0) then begin
- WIDGET_CONTROL, MenuButtons(i), SENSITIVE=0
- endif
-
- if (flag and 1) then begin
- level = level + 1
- parent(level) = MenuButtons(i)
- names(level) = txt + '|'
- lflags(level) = (flag AND 2) ne 0
- endif else if ((flag and 2) NE 0) then begin
-
- ; Pop the previous levels.
- ;
- while (lflags(level)) do level = level-1
-
- ; Pop this level.
- ;
- level = level - 1
- endif
- endfor
-
- end ; of DemoMenuCreate
-
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Create the reconstruction filter : none
- ;
- function None, $
- x, $ ; IN: x coordinates
- pointSpacing ; IN: point spacing
-
- RETURN, [1.0]
- end ; of None
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Create the reconstruction filter : RAMLAK
- ;
- function Ramlak, $
- x, $ ; IN: x coordinates
- pointSpacing ; IN: point spacing
-
- zero = where(x EQ 0.0, count)
- q = x
- if (count NE 0) then q(zero) = .01
- y = -SIN(!pi*x/2)^2 / (!pi^2 * q^2 * pointSpacing)
- if (count NE 0) then y(zero) = 1./(4.*pointSpacing)
- RETURN, y
- end
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Create the reconstruction filter : Shepp_Logan
- ;
- function Shepp_Logan,$
- x, $ ; IN: x coordinates
- pointSpacing ; IN: point spacing
-
- pointSpacing = !pi^2 * pointSpacing * (1.-4.*x^2)
- zeros = where(abs(pointSpacing) LE 1.0e-6, count)
- if (count ne 0) then pointSpacing(zeros) = .001
- RETURN, 2./pointSpacing
- end
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Create the reconstruction filter : Lp_Cosine
- ;
- function Lp_Cosine, $
- x, $ ; IN: x coordinates
- pointSpacing ; IN: point spacing
-
- RETURN, 0.5 * (ramlak(x-.5, pointSpacing) + $
- ramlak(x+.5, pointSpacing))
- end
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Create the reconstruction filter : Gen_Hamming
- ;
- function Gen_Hamming, $
- x, $ ; IN: x coordinates
- pointSpacing, $ ; IN: point spacing
- alpha ; IN: alpha attenuation factor
-
- if (N_ELEMENTS(alpha) LE 0) then alpha = 0.5
- RETURN, alpha * ramlak(x, pointSpacing) + ((1.-alpha)/2) * $
- (ramlak(x-1, pointSpacing) + ramlak(x+1, pointSpacing))
- end
-
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Fill an elliptical area within an array with a given value.
- ;
- pro DEllipse, $
- regionArray, $ ; IN/OUT: array
- x0, $ ; IN: x loaction of the ellipse center
- y0, $ ; IN: y loaction of the ellipse center
- rx, $ ; IN: radius in x of the ellipse
- ry, $ ; IN: radius in y of the ellipse
- theta, $ ; IN: angle of the major axis
- value ; IN: value given to a point within the ellipse.
-
- s = SIZE(regionArray)
- nx = s(1)
- ny = s(2)
- n = 64 ;# of points around ellipse
- x = FINDGEN(n) * (2 * !pi / n)
- y = ry * SIN(x)
- x = rx * COS(x)
- t = theta * !dtor
- yp = (COS(t) * y + SIN(t) * x + y0) * (ny/2) + (ny/2)
- x = (COS(t) * x - SIN(t) * y + x0) * (nx/2) + (nx/2)
- regionArray(POLYFILLV(x, yp, nx, ny)) = value
- end
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Fill a circular area within an array with a given value.
- ;
- pro DCircle, $
- regionArray, $ ; IN/OUT: array
- r, $ ; IN: radius of the circle
- x0, $ ; IN: x location of the circle center
- y0, $ ; IN: y location of the circle center
- value ; IN; value
-
- s = SIZE(regionArray)
- n = 100
- x = FINDGEN(n) * (2 * !pi / n)
- y = r * SIN(x) + y0
- x = r * COS(x) + x0
- regionArray(POLYFILLV(x, y, s(1), s(2))) = value
- end
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Draw a polygon within an array
- ;
- pro DPoly, $
- regionArray, $ ; IN/OUT: array
- x, $ ; IN: x location of the polygon center
- y, $ ; IN: y location of the polygon center
- value ; IN: value
-
- s =size(regionArray)
- regionArray(POLYFILLV(x, y, s(1), s(2))) = value
- end
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Create the phantom (image) objects.
- ;
- function Make_Phantom, $
- imageSize, $ ; IN: Image size
- object ; IN: Object identifier, see list below
-
- case object of
-
- ; Create squares.
- ;
- 0: Begin
- rslt = FLTARR(imageSize, imageSize)
- rslt(imageSize/2, imageSize/2) = 1.0
- endcase
-
- ; Create ellipses.
- ;
- 1: begin
- rslt = FLTARR(imageSize, imageSize)
- dellipse, rslt, 0, 0, .69, .92, 0, 1.0 ;a
- dellipse, rslt, 0,-.0184, .6624, .8740, 0, 0.34 ;b
- dellipse, rslt, .22, 0., .11, .31, -18, 0. ;c
- dellipse, rslt, -.22, 0., .16, .41, 18, 0. ;d
- dellipse, rslt, 0., .35, .21, .25, 0., 0.67 ;e
- dellipse, rslt, 0., .35, .05, .05, 0., 0.57 ;Extra
- dellipse, rslt, 0., .1, .046, .046, 0., .67 ;f
- dellipse, rslt, 0., -.1, .046, .046, 0., .67 ;g
- dellipse, rslt, -.08, -.605, .046, .023, 0., .67 ;h
- dellipse, rslt, 0., -.606, .023, .023, 0., .67 ;i
- dellipse, rslt, 0.06, -.605, .023, .046, 0., .67 ;j
- dellipse, rslt, -.49, -.470, .050, .050, 0., 1.00 ;j
- endcase
-
- ; Create circles.
- ;
- 2: begin
- rslt = FLTARR(imageSize, imageSize)
- dcircle, rslt, imageSize/4, .5*imageSize, .5*imageSize, .5
- dcircle, rslt, imageSize/12, .6*imageSize, .6*imageSize, 0.0
- dcircle, rslt, imageSize/12, .38*imageSize, .37*imageSize, 1.
- dcircle, rslt, imageSize/12, .7*imageSize, .2*imageSize, 1.
- dcircle, rslt, 3, .28*imageSize, .60*imageSize, 0.0
- dpoly, rslt, imageSize*[.55, .60, .65], imageSize * [.37, .45, .37], 0
- rslt(imageSize*[.117, .195], imageSize*[.781, .855]) = .5
- endcase
-
- ; Create polygons.
- ;
- 3: begin
- rslt = FLTARR(imageSize, imageSize)
- x = imageSize/3
- y = 2 * imageSize / 3
- s = 2
- rslt(x-s:x+s, y-s:y+s) = 1.0
- rslt(y,x)=1.0
- s = 4
- rslt(x-s:x+s, x-s:x+s) = 0.5
- dpoly, rslt, [y-s, y+s, y], [y-s, y-s, y+s], 1.0
- endcase
-
- ; Download the CT scan image (slice), resize, and scale it.
- ;
- 4: begin
- OPENR,unit, filepath('ctscan.dat', $
- SUBDIR=['examples','data']), $
- /GET_LUN
- imageArray=BYTARR(256,256)
- READU, unit, imageArray
- CLOSE, unit
- FREE_LUN, unit
- imageArray = (imageArray < 200b)/ 200. ;Normalize
- if (imageSize LT 256) then rslt = CONGRID(imageArray, $
- imageSize, imageSize, /INTERP) $
- else if (imageSize GT 256) then begin
- rslt = FLTARR(imageSize, imageSize)
- rslt((imageSize-256)/2, (imageSize-256)/2) = imageArray ;Insert
- endif else rslt = imageArray ;already correct size.
- endcase
- endcase
-
- RETURN, rslt
-
- end ; of Make_Phantom
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Reconstruct (recompute) the image and display it.
- ;
- pro Reconstruct_It
-
- COMMON tomodemo_com, base, imageSize, $
- interp, nangles, kernelSize, filter, $
- draw, window, labels, lnames, top, shiftValue, pointSpacing, $
- filters, obj_button, object, zoomFactor, maxa, mina, ocolors, $
- MenuItems, MenuButtons, wReconstructButton, wAngleButton, wFilterButton, $
- wKernelButton, wInterpButton, drawErrorID, errorBase, $
- Bar_Base, wSubBase, $ ; base ID to desensitize when computing.
- ReconstructFlag, $ ; 0 = Go button has not been pushed yet
- originalImage, reconImage, sinogramImage, $ ; Images
- sText, wText, $ ; widget ids and text structure for tips.
- pimage, pimagestr, nimages, comp_view ;history
-
- ; Desensitize the bases during reconstruction.
- ;
- WIDGET_CONTROL, Bar_base, SENSITIVE=0
- WIDGET_CONTROL, wSubBase, SENSITIVE=0
-
- ; Initialize the windows labels (names).
- ;
- lnames2 = ['Original (Click window)', $
- 'Reconstruction ', $
- 'Sinogram (Click window) ', $
- 'Error (Click window) ']
- ; Create the widgets starttin g with the top level base.
- for i=0,3 do begin
- ; WIDGET_CONTROL, labels(i), SET_VALUE=lnames(i)
- WIDGET_CONTROL, labels(i), SET_VALUE=lnames2(i)
- endfor
-
- ; Redraw the original image by showing the angle line.
- ;
- WSET, window(0)
- TV, BYTSCL(originalImage, TOP=top, MAX=maxa, MIN=mina)
-
- ; Draw a circle.
- ;
- x = FINDGEN(60) * (!pi * 2 / 59)
- na2 = imageSize/2.
- PLOTS, na2 * COS(x) +na2, na2*SIN(x) +na2, /DEVICE, COLOR=top+1
- comp_view = 0
-
- ; Sort indices by i,j, where k = nangles * j / 2^i
- ;
- ind = INDGEN(nangles) ;Indices of projections
- if (0) then begin
- j=1
- mm = CEIL(ALOG(nangles) / ALOG(2.0))
-
- for k=1, mm do for i=1,2L^k-1,2 do begin
- kk = i*2L^mm/2L^k
- if (kk LT nangles) then begin
- ind(j) = kk
- j = j+1
- endif
- endfor
- endif
-
- ; Compute the filter function.
- ;
- x = FINDGEN(2*kernelSize+1)-kernelSize
- ker = call_function(filters(filter), x, 1.0) ; Point spacing is 1.0
- ker = ker / (TOTAL(ker) * nangles * SQRT(imageSize)) ; Normalize it
-
- np = FIX(SQRT(2.) * imageSize/pointSpacing + kernelSize + 4)
- sinogramImage = FLTARR(np, nangles)
- convolvedSinogramImage = FLTARR(np, nangles)
- t0 = 0.
- shiftValue = -(np-imageSize)/2
- WIDGET_CONTROL, labels(3), SET_VALUE='Convolved Sinogram '
- chunk = nangles/32 > 1
- first = 1
-
- ; Draw the image section for every angle.
- ;
- for i=0,nangles-1 do begin
- WSET, window(0)
- DEVICE, SET_GRAPH=6
-
- ; Erase the old line.
- ;
- if (i NE 0) then $
- PLOTS, [na2+x, na2-x], [na2+y, na2-y], /DEVICE, COLOR=top/2
-
- t = i * !pi / nangles
- t1 = SYSTIME(1)
- x = SIN(t) * na2
- y = COS(t) * na2
-
- ; Draw the line except for last one.
- ;
- if (i NE (nangles-1)) then $
- PLOTS, [na2+x, na2-x], [na2+y, na2-y], /DEVICE, COLOR=top/2
-
- ; Compute the various images.
- ;
- device, SET_GRAPH=3
- p = FLTARR(np)
- riemann, p, originalImage, t, BILINEAR=interp EQ 1, $
- CUBIC=interp EQ 2, D=pointSpacing
- sinogramImage(0,i) = p
- convolvedSinogramImage(0,i) = CONVOL(p, ker)
- t0 = t0 + SYSTIME(1) - t1
-
- if ((i LT imageSize) AND ((i mod chunk) EQ (chunk-1)) ) then begin
- WSET, window(2)
- if (first) then erase
- i0 = i - chunk+1
- TV, BYTSCL(shift(sinogramImage(*,i0:i), shiftValue), TOP=top), 0, i0
- WSET, window(3)
- if (first) then erase
- first = 0
- TV, BYTSCL(SHIFT(convolvedSinogramImage(*,i0:i), shiftValue), $
- TOP=top), 0, i0
- endif
-
- str = "View "+STRTRIM(i,2)
- sText.text[6] = str
- textChange = ['views','void1']
- putTips, sText, wText[1], $
- textChange, [1,2]
-
- endfor
-
- str = 'execution time : ' + string(t0, FORMAT='(f7.2)') + ' sec.'
- sText.text[7] = str
- textChange = ['xtime']
- putTips, sText, wText[1], $
- textChange, [2]
-
- ; Display the sinogram.
- ;
- WSET, window(2)
- sinogramImage = BYTSCL((SHIFT(sinogramImage, $
- shiftVAlue))(0:imageSize-1,*), TOP=top)
- zoomFactor = 1 ;Zoom factor
- TV, sinogramImage
-
- ; Draw the convolved sinogram image.
- ;
- WSET, window(3)
- TV, BYTSCL(shift(convolvedSinogramImage, shiftValue), TOP=top)
-
- ; Compute and draw the reconstructed image.
- ;
- reconImage = FLTARR(imageSize, imageSize)
-
- step = 1 > (nangles/32)
- t0 = 0.
-
- WSET, window(1)
- erase
- for i=0,nangles-1 do begin
- j = ind(i)
- t = j * !pi / nangles ;The angle
- t1 = SYSTIME(1)
-
- riemann, convolvedSinogramImage(*,j), $
- reconImage, t, /BACK, D=pointSpacing, $
- BILINEAR=interp EQ 1, CUBIC=interp EQ 2
-
- t0 = t0 + SYSTIME(1)-t1
- if (i MOD step EQ 0) then begin
-
- if 0 then begin ;NO OPPED
- WSET, window(0)
- DEVICE, SET_GRAPH=6
- if i ne 0 then $
- PLOTS, [na2+x, na2-x], [na2+y, na2-y], /DEVICE, COLOR=top/2
- x = -SIN(t) * na2
- y = COS(t) * na2
- if i ne (nangles-1) then $ ;Erase old
- PLOTS, [na2+x, na2-x], [na2+y, na2-y], /DEVICE, COLOR=top/2
- DEVICE, SET_GRAPH=3
- WSET, window(1)
- endif ;NO OPPED
-
- TV, BYTSCL(reconImage, TOP=top)
- endif
-
- str = "View "+STRTRIM(i,2)
- sText.text[6] = str
- textChange = ['views']
- putTips, sText, wText[1], $
- textChange, [1]
-
- endfor
-
- scaledReconImage = BYTSCL(reconImage, TOP=top)
- tv, scaledReconImage
-
- meana = TOTAL(originalImage)/N_ELEMENTS(originalImage)
- meanb = TOTAL(reconImage)/N_ELEMENTS(reconImage)
- ss = SQRT(TOTAL((originalImage-meana)^2) / TOTAL((reconImage-meanb)^2))
- reconImage = (reconImage-meanb) * ss + meana
- errorImage = abs(originalImage-reconImage)
- errtot = LONG(TOTAL(errorImage^2))
-
- str1 = 'execution time : ' + string(t0, FORMAT='(f7.2)') + ' sec.'
- str2 = "error norm : "+ string(errtot, FORMAT="(i8)")
- sText.text[7] = str1
- sText.text[8] = str2
- textChange = ['error','xtime']
- putTips, sText, wText[1], $
- textChange, [1,2]
-
- ; Draw the error image.
- ;
- WSET, window(3)
- WIDGET_CONTROL, labels(3), $
- SET_VALUE='Error (Click window)'
- TV, BYTSCL(errorImage, MAX=maxa, MIN=0, TOP=top)
-
- if (nimages EQ 0) then begin ;Save history
- pimage = BYTARR(imageSize, imageSize, 4, /NOZERO)
- pimagestr = STRARR(4)
- endif
-
- pimage(0, 0, nimages MOD 4) = scaledReconImage
- pimagestr(nimages MOD 4) = STRING(filters(filter), STRTRIM(nangles,2), $
- ([ 'NN', 'In','Cu'])(interp), STRTRIM(kernelSize*2+1,2), STRTRIM(errtot,2),$
- FORMAT="(a4, ',V',a, ',',a, ',K',a,',E',a)")
- nimages = nimages+1
-
- ; Resensitze the bases when reconstruction is done.
- ;
- WIDGET_CONTROL, Bar_base, SENSITIVE=1
- WIDGET_CONTROL, wSubBase, SENSITIVE=1
-
- end ; reconstruct_it
-
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Redraw the phanton, sinogram, etc.
- ;
- pro Recon_Redraw, $
- force ; IN: flag indicator
-
- common tomodemo_com
-
- ; Do nothing (return) under certain conditions.
- ;
- if ( (NOT comp_view) AND (force eq 0)) then RETURN
-
- ; Reset the original window titles.
- ;
- for i= 0,3 do WIDGET_CONTROL, labels(i), SET_VALUE=lnames(i)
-
- WSET, window(0)
- TV, BYTSCL(originalImage, TOP=top, MIN=mina, MAX=maxa)
-
- if (ReconstructFlag EQ 0) then begin
- WSET, window(2)
- erase
- WSET, window(1)
- erase
- WSET, window(3)
- erase
- endif
-
- if (N_ELEMENTS(sinogramImage) GT 1) then begin
- WSET, window(2)
- erase
-
- if (force EQ 1) then begin
- TV, REBIN(sinogramImage, imageSize, zoomFactor*nangles)
- endif
-
- WSET, window(1)
- erase
-
- if (force EQ 1) then begin
- TV, BYTSCL(reconImage, TOP=top)
- endif
-
- WSET, window(3)
- erase
-
- if (force EQ 1) then begin
- TV, BYTSCL(ABS(originalImage-reconImage), MAX=maxa, MIN=0, TOP=top)
- endif
- endif
-
- textChange = ['error','xtime']
- putTips, sText, wText[1], $
- textChange, [1,2]
-
- comp_view = 0
- end ; of recon_redraw
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Main event handler
- ;
- pro Rec_Demo_Event, $
- sEvent ; IN: event structure
-
- common tomodemo_com
-
- ; Quit the application using the close box.
- ;
- if (TAG_NAMES(sEvent, /STRUCTURE_NAME) EQ $
- 'WIDGET_KILL_REQUEST') then begin
- WIDGET_CONTROL, sEvent.top, /DESTROY
- RETURN
- endif
-
- forward_function DemoMenuChoice
-
- WIDGET_CONTROL, sEvent.top, /HOURGLASS
-
- ; Branch accordingly to the event ID
- ; (which widget created that event)
- ;
- case sEvent.id of
-
- ; Reconstruct the image.
- ;
- wReconstructButton: goto, do_reconstruction
-
- ; Button press within the original view area.
- ; Draw the cursor location with a small square on the original
- ; image and draw the corresponding curve on the sinogram.
- ;
- draw(0): begin
-
- ; Return if the button was not the left mouse button or
- ; the composite view does not exist or the
- ; sinogram image does not exist.
- ;
- if ((sEvent.press NE 0) OR comp_view) then RETURN
- if (N_ELEMENTS(sinogramImage) LE 1) then RETURN
-
- ocolors(0) = (ocolors(0)+1) mod 4
- c = ocolors(0) + top + 1
- WSET, window(0)
- PLOTS, sEvent.x, sEvent.y, COLOR=c, /DEVICE, PSYM=6
- x = (sEvent.x - imageSize/2)
- y = (sEvent.y - imageSize/2)
- t = FINDGEN(nangles+1)* (!pi / nangles) ;Rotate by theta
- t = (COS(t) * x + SIN(t) * y) + (imageSize/2)
-
- ; Draw the line on the sinogram image.
- ;
- WSET, window(2)
- i = (imageSize /(2* nangles)) > 1
-
- if (i NE zoomFactor) and (N_ELEMENTS(sinogramImage) GT 1) then begin
- zoomFactor = i
- TV, REBIN(sinogramImage, imageSize, zoomFactor * nangles)
- endif
-
- PLOTS, t, zoomFactor*FINDGEN(nangles), COLOR=c, /DEVICE
-
- empty
-
- endcase ; of draw(0)
-
- ; Show the error on the horizontal line where the
- ; left mouse button was pressed on the error view.
- ; Open an error view window if it does not already exist.
- ;
- draw(3): begin
-
- ; Return if the button was not the left mouse button or
- ; the composite view does not exist or the
- ; reconstructed image does not exist or the
- ; error or sinogram images do not exist.
- ;
- if ((sEvent.press Ne 0) OR comp_view) then RETURN
- if (N_ELEMENTS(reconImage) LE 1) then RETURN
- if (N_ELEMENTS(sinogramImage) LE 1) then RETURN
-
- charscale = 8.0/!d.X_CH_SIZE
-
- ; Check if the error plot already exist, if so return
- ;
- errorStatus = WIDGET_INFO(drawErrorID, /VALID_ID)
- if (errorStatus EQ 0) then begin
-
- ; Create a new window that displays the error plot.
- ;
- errorBase = WIDGET_BASE(TLB_FRAME_ATTR=1, $
- TITLE ='Error', $
- XOFFSET=75, YOFFSET=75, $
- GROUP_LEADER=base)
-
- drawErrorID = WIDGET_DRAW(errorBase, $
- SCR_XSIZE=250, SCR_YSIZE=250)
-
- WIDGET_CONTROL, errorBase, /REALIZE
- endif
-
- WIDGET_CONTROL, drawErrorID, GET_VALUE=errorWindow
-
- WSET, errorWindow
- y = sEvent.y
- PLOT, originalImage(*,y), YRANGE=[mina, maxa], XSTYLE=3, YSTYLE=2, $
- TITLE='Row '+STRTRIM(y,2), XMARGIN=[4,1], YMARGIN=[2,2], $
- CHARSIZE = 1.0*charscale
-
- OPLOT, reconImage(*,y), COLOR=top+2, PSYM=3
-
- empty
-
- endcase ; of draw(3)
-
- ; When the left mouse button is pressed within the
- ; sinogram view area, draw a small square to mark the
- ; cursor location and draw the corresponding curve on the
- ; original image.
- ;
- draw(2): begin
-
- ; Return if the button was not the left mouse button or
- ; the composite view does not exist or the
- ; sinogram image does not exist.
- ;
- if ((sEvent.press NE 0) OR comp_view) then RETURN
- if N_ELEMENTS(sinogramImage) le 1 then return
- if N_ELEMENTS(zoomFactor) le 0 then zoomFactor = 1 ;Within image?
-
- if (sEvent.y GE nangles*zoomFactor) then RETURN
-
- ocolors(2) = (ocolors(2) + 1) mod 4
- c = ocolors(2) + top + 1
- na2 = imageSize/2
- WSET, window(2)
- PLOTS, sEvent.x, sEvent.y, COLOR=c, /DEVICE, PSYM=6
-
- s = sEvent.x - na2 ; Radial distance
- t = sEvent.y/zoomFactor * !pi / nangles ; Angle in radians
-
- ; Display the angle and the x distance from center
- ; informations.
- ;
- str1 = 'Theta :' + string(t*180./!pi, FORMAT="(f5.1)")
- str2 = 'X from center :' + string(s, FORMAT="(f6.1)")
- sText.text[10] = str1
- sText.text[11] = str2
- textChange = ['theta', 'cente']
- putTips, sText, wText[1], $
- textChange, [1,2]
-
- x = [-na2, 0, na2]
- st = SIN(t)
- ct = COS(t)
-
- if (t NE 0.0) then begin
- y = (s - x * ct) / st + na2
- endif else begin
- y=[0, imageSize-1]
- x = [sEvent.x, sEvent.x]-na2
- endelse
-
- ; Draw the corresponding curve on the original image.
- ;
- WSET, window(0)
- PLOTS, x+na2, y, /DEVICE, COLOR=c ;Show projection
- py = sinogramImage(*, sEvent.y/zoomFactor) ;Draw the profile, rotated....
- py = py * (imageSize/(5.*max(py)))
- px = FINDGEN(imageSize) - na2
- PLOTS, ct * px - st * py + na2, $
- st * px + ct * py + na2 , /DEVICE, COLOR=c
- empty
- endcase ; of draw(2)
-
- ; Do nothing (return) when the mouse button is pressed
- ; on the reconstructed image.
- ;
- draw(1): RETURN
-
- ; Set the number of angles.
- ;
- wAngleButton: nangles = 2^(sEvent.index+2)
-
- ; Set the type of filter.
- ;
- wFilterButton: filter = sEvent.index
-
- ; Set the kernel size.
- ;
- wKernelButton: kernelSize = 2^sEvent.index
-
- ; Set the interpolation method.
- ;
- wInterpButton: interp = sEvent.index
-
- ; Any other event must be a menu button.
- ;
- else: begin
-
- ; Get the user value of the button.
- ;
- WIDGET_CONTROL, sEvent.id, GET_UVALUE=uv
-
- uv1 = STR_SEP(uv, "|")
-
- ; Branch to the appropriate button event.
- ;
- case uv1(1) of
- 'File': case uv1(2) of
-
- ; Quit this application.
- ;
- 'Quit': begin
- originalImage = 0 ;Remove images & clean up.
- reconImage=0
- sinogramImage = 0
- pimage=0
- pimagestr=0
- WIDGET_CONTROL, sEvent.top, /DESTROY
- RETURN
- endcase
-
- ; Draw the newly selected object and destroy the
- ; previous one and its associated images.
- ;
- 'Object': begin
- WSET, window(0)
- object = DemoMenuChoice(uv, MenuItems, MenuButtons) + 1
- originalImage = make_phantom(imageSize, object)
- maxa = max(originalImage, MIN=mina)
- sinogramImage = 0
- ReconstructFlag = 0
- recon_redraw, 1
-
- ; Destroy the error window.
- ;
- if ( widget_info(errorbase, /VALID_ID)) then begin
- WIDGET_CONTROL, errorBase, /DESTROY
- endif
- nimages = 0
-
- ; Set the execution time, error norm, theta, and
- ; center text fields to nothing.
- ;
- sText.text[7] = ' '
- sText.text[8] = ' '
- sText.text[10] = ' '
- sText.text[11] = ' '
-
- ; Reset the tips text.
- ;
- textChange = ['click', 'numer']
- putTips, sText, wText[1], $
- textChange, [1,2]
-
- endcase
-
- 'Reconstruct': begin
- do_reconstruction:
- WIDGET_CONTROL, base, /HOURGLASS
- reconstruct_it
- ReconstructFlag = 1
-
- ; Destroy the error window.
- ;
- if ( widget_info(errorbase, /VALID_ID)) then begin
- WIDGET_CONTROL, errorBase, /DESTROY
- endif
- endcase
-
- endcase ; of uv1(2) or File
-
- ; Display the information text file.
- ;
- 'About': begin
- if (Xregistered('Xdisplayfile') NE 0) then RETURN
-
- Xdisplayfile, filepath("reconstr.txt", $
- SUBDIR=['examples','demo','demotext']), $
- TITLE="About the Reconstruction Demo", $
- GROUP=sEvent.top, DONE_BUTTON='Done', $
- WIDTH=55, HEIGHT=14
- endcase
-
- ; Load a new color table.
- ;
- 'Edit': begin
- if (Xregistered('XLoadct') NE 0) then RETURN
- XLOADCT, NCOLORS=top, GROUP=sEvent.top
- endcase
-
-
- 'View': case uv1(2) of
-
- ; Redraw all the views.
- ;
- 'Redraw': begin
- recon_redraw, 1
- endcase
-
- ; Compare the previous reconstructed image with the
- ; current one for the same object.
- ;
- 'Compare': begin
- if (nimages le 1) then return
- n = nimages < 4
- for i=nimages-n, nimages-1 do begin
- j = i - nimages+n
- WSET, window(j)
- WIDGET_CONTROL, labels(j), SET_VALUE=pimagestr(i mod 4)
- TV, pimage(*,*,i mod 4)
- endfor
- comp_view = 1
- endcase
-
- endcase ; of View
-
- endcase ; of uv1 (menu bar buttons)
-
- endcase ; of else
-
- endcase ; of sEvent.id
-
- end ; of event handler
-
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Cleanup procedure.
- ;
- pro Reconstr_Cleanup, $
- tlb ; IN: top level base identifier
-
- ; Get the color table saved in the window's user value.
- ;
- WIDGET_CONTROL, tlb, GET_UVALUE=info, /NO_COPY
-
- ; Restore the previous color table.
- ;
- TVLCT, info.colorTable
-
- ; Map the group leader base if it exists.
- ;
- if (WIDGET_INFO(info.groupBase, /VALID_ID)) then $
- WIDGET_CONTROL, info.groupBase, /MAP
-
- end ; Of Reconstr_Cleanup
-
- ;--------------------------------------------------------------------
- ;
- ; PURPOSE Cleanup procedure.
- ;
- pro D_Reconstr, $
- ImageSizeIn, $ ; IN: (opt) image size vector
- GROUP=group, $ ; IN: (opt) group identifier
- APPTLB = appTLB ; OUT: (opt) TLB of this application
-
- common tomodemo_com
-
- ; Check the validity of the group identifier.
- ;
- ngroup = N_ELEMENTS(group)
- if (ngroup NE 0) then begin
- check = WIDGET_INFO(group, /VALID_ID)
- if (check NE 1) then begin
- print,'Error, the group identifier is not valid'
- print, 'Return to the main application'
- RETURN
- endif
- groupBase = group
- endif else groupBase = 0L
-
- ; Save the current color table.
- ;
- TVLCT, savedR, savedG, savedB, /GET
- colorTable = [[savedR], [savedG], [savedB]]
-
- ; Initialize sizes.
- ;
- if (N_ELEMENTS(ImageSizeIn) LE 0) then begin
- DEVICE, GET_SCREEN = x
- if (x(0) GT 1024) then ImageSize=256 $
- else if x(0) ge 800 then ImageSize=192 $
- else ImageSize = 128
- endif else begin
- imageSize = imageSizeIn
- endelse
-
- ; Get the tips.
- ;
- sText = getTips(filepath('reconstr.tip', $
- SUBDIR=['examples','demo', 'demotext']) )
-
- ; Create the starting up message.
- ;
- if (ngroup EQ 0) then begin
- drawbase = startmes()
- endif else begin
- drawbase = startmes(GROUP=group)
- endelse
-
-
- ; Initialize working parameters and arrays.
- ;
- imageSize = ImageSize ; Default values
- filter = 2 ; Shepp Logan
- interp = 1 ; Linear
- nangles = 32
- kernelSize = 8
- pointSpacing = 1. ; Point spacing
- zoomFactor = 1
- object = 4
- n = 2 * imageSize
- draw = LONARR(4)
- window = LONARR(4)
- ocolors = INTARR(4) ; Overlay colors
- pimage=0
- nimages=0
- comp_view = 0 ; Showing normal view
- ReconstructFlag = 0
-
- ; Initialize the filter options, None must be 1st filter....
- ;
- filters= ['None', 'RamLak', 'Shepp_Logan', $
- 'LP_Cosine', 'Gen_Hamming']
-
- ; Initialize the windows labels (names).
- ;
- lnames = ['Original ', $
- 'Reconstruction ', $
- 'Sinogram ', $
- 'Error ']
- ; Create the widgets starttin g with the top level base.
- ;
- if (N_ELEMENTS(group) EQ 0) then begin
- base = WIDGET_BASE(TITLE='Reconstruction Demo', $
- MAP=0, $
- /TLB_KILL_REQUEST_EVENTS, $
- MBAR=bar_base, TLB_FRAME_ATTR=1, /COLUMN)
- endif else begin
- base = WIDGET_BASE(TITLE='Reconstruction Demo', $
- GROUP_LEADER=group, $
- MAP=0, $
- /TLB_KILL_REQUEST_EVENTS, $
- MBAR=bar_base, TLB_FRAME_ATTR=1, /COLUMN)
- endelse
-
- MenuItems = [ $
- '1File', $
- '1Object', '0Shepp-Logan Phantom', '0Circles', '0Squares','6CT Slice', $
- '0Reconstruct', '2Quit', $
- '1Edit', '2Color Palette', $
- '1View', '0Redraw', '2Compare', $
- '1About', '2About Reconstruction' ]
-
- ; Create the menu bar and all its buttons.
- ;
- DemoMenuCreate, MenuItems, MenuButtons, Bar_base
-
- ; Create a sub base that has the left and right bases.
- ;
- wSubBase = WIDGET_BASE(base, COLUMN=2)
-
- ; Create the base that has options and functionality buttons.
- ;
- wLeftBase = WIDGET_BASE(wSubBase, /COLUMN, /BASE_ALIGN_CENTER)
-
- ; Create the angel base.
- ;
- wViewBase = WIDGET_BASE(wLeftBase, /COLUMN)
-
- ; Create the angel base.
- ;
- wAngleBase = WIDGET_BASE(wViewBase, $
- /BASE_ALIGN_LEFT, /COLUMN)
-
- ; Create the number of angles droplist.
- ;
- wAnglesLabel = WIDGET_LABEL(wAngleBase, $
- VALUE='Number of angles :')
-
- ; Droplist to select the nuimber of angles.
- ;
- wAngleButton = WIDGET_DROPLIST(wAngleBase, $
- VALUE= ['4','8','16','32','64','128','256'])
-
- ; Create the filter base.
- ;
- wFilterBase = WIDGET_BASE(wLeftBase, /COLUMN, /FRAME)
-
- ; Create the number of angles droplist.
- ;
- wFilterLabel = WIDGET_LABEL(wFilterBase, $
- /ALIGN_CENTER, $
- VALUE='Filter')
-
- ; Create the sub base of filter base.
- ;
- wSubFilterBase = WIDGET_BASE(wFilterBase, $
- /BASE_ALIGN_LEFT, /COLUMN)
-
- ; Create the filter type label.
- ;
- wTypeLabel = WIDGET_LABEL(wSubFilterBase, $
- VALUE='Type :')
-
- ; Droplist to select the filter type.
- ;
- wFilterButton = WIDGET_DROPLIST(wSubFilterBase, $
- VALUE=filters)
-
- ; Create the kernel size label.
- ;
- wTypeLabel = WIDGET_LABEL(wSubFilterBase, $
- VALUE='Kernel size :')
-
- ; Droplist to select the kernel size.
- ;
- wKernelButton = WIDGET_DROPLIST(wSubFilterBase, $
- VALUE=['3','5','9','17','33','65'])
-
- ; Create the interpolation label.
- ;
- wInterpolationLabel = WIDGET_LABEL(wSubFilterBase, $
- VALUE='Interpolation :')
-
- ; Droplist to select the interpolation method.
- ;
- wInterpButton = WIDGET_DROPLIST(wSubFilterBase, $
- VALUE=['Nearest Neighbor','Linear','Cubic'])
-
-
- ; Create the reconstruct button.
- ;
- wReconstructButton = WIDGET_BUTTON(wLeftBase, $
- VALUE='Reconstruct', /NO_RELEASE)
-
- ; Create the base that has options and functionality buttons.
- ;
- wRightBase = WIDGET_BASE(wSubBase, /COLUMN)
-
- wRow1Base = WIDGET_BASE(wRightBase, /ROW) ;2 x 2 bases
- temp = LONARR(4)
- labels = LONARR(4)
-
- for i=0,1 do temp(i) = WIDGET_BASE(wRow1Base, /COLUMN)
-
- wRow2Base = WIDGET_BASE(wRightBase, /ROW)
-
- for i=2,3 do temp(i) = WIDGET_BASE(wRow2Base, /COLUMN)
-
- for i=0,3 do begin
- labels(i) = WIDGET_LABEL(temp(i), $
- /ALIGN_LEFT, $
- /DYNAMIC_RESIZE, $
- VALUE=STRING(lnames(i), FORMAT='(a34)'))
-
- draw(i) = WIDGET_DRAW(temp(i), /BUTTON, $
- RETAIN=2, $
- XSIZE=imageSize, YSIZE=imageSize)
- endfor
-
- for i=0,3 do WIDGET_CONTROL, labels(i), SET_VALUE=lnames(i)
-
- ; Create the status line label.
- ;
- wStatusBase = WIDGET_BASE(base, MAP=0, /ROW)
-
- nWidgets = 2
- wText = LONARR(nWidgets)
- widTips, wStatusBase, sText.text, XSIZE=36, $
- YSIZE=3, NWIDGETS=nWidgets, wText
-
-
-
- ; Realize the widget hierarchy.
- ;
- WIDGET_CONTROL, base, /REALIZE
-
- WIDGET_CONTROL, wAngleButton, SET_DROPLIST_SELECT=3
- WIDGET_CONTROL, wFilterButton, SET_DROPLIST_SELECT=filter
- WIDGET_CONTROL, wKernelButton, SET_DROPLIST_SELECT=3
- WIDGET_CONTROL, wInterpButton, SET_DROPLIST_SELECT=1
-
- ; Make the angle and filter base the same x dimension.
- ;
- szAngle = WIDGET_INFO(wAngleBase, /GEOMETRY)
- szFilter = WIDGET_INFO(wSubFilterBase, /GEOMETRY)
- xAngle = szAngle.scr_xsize + ( 2* szAngle.margin)
- xFilter = szFilter.scr_xsize + ( 2* szFilter.margin)
- if (xAngle LT xFilter) then begin
- WIDGET_CONTROL, wAngleBase, SCR_XSIZE=xFilter
- endif else begin
- WIDGET_CONTROL, wFilterBase, SCR_XSIZE=xAngle
- endelse
-
- ; Returns the top level base to the APPTLB keyword.
- ;
- appTLB = base
-
- ; Size the tips widgets.
- ;
- sizeTips, base, wText, wStatusBase
-
- ; Get the windows (drawing areas) identifiers.
- ;
- for i=0,3 do begin
- WIDGET_CONTROL, draw(i), GET_VALUE=j
- window(i) = j
- endfor
- top = !D.TABLE_SIZE-6
-
- ; Load the grey scale colr table.
- ;
- LOADCT, 0, /SILENT, NCOLORS=top+1
-
- ; Allocate working colors : red, green, yellow, blue, white.
- ;
- TVLCT, [255,0,255,0,255],[0,255,255,0,255],[0,0,0,255,255], top+1
-
- ; Destroy the starting up window.
- ;
- WIDGET_CONTROL, drawbase, /DESTROY
-
- ; Map the top level base.
- ;
- WIDGET_CONTROL, base, MAP=1
-
- ; Load the CT scan image as default.
- ;
- originalImage = make_phantom(imageSize, object)
- maxa = max(originalImage, MIN=mina)
- WSET, window(0)
- TV, BYTSCL(originalImage, TOP=top, MAX=maxa, MIN=mina)
-
- ; Assign an initial value to drawErrorID.
- ;
- drawErrorID = -1L
- errorBase = -1L
-
- ; Create the info structure.
- ;
- info = { $
- ColorTable: colorTable, $
- groupBase: groupBase $ ; Base of Group Leader
- }
-
- ; Assign the info structure into the user value of the top level base.
- ;
- WIDGET_CONTROL, base, SET_UVALUE=info, /NO_COPY
-
- XMANAGER, 'd_reconstr', base, CLEANUP='Reconstr_Cleanup', $
- Event_Handler='rec_demo_event', /NO_BLOCK
- end
-