Rapid geometry interrogation for a uniform volume element-based Monte Carlo particle transport simulation by Michael William Frandsen A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science Montana State University © Copyright by Michael William Frandsen (1998) Abstract: New particle tracking methods to speed Boron Neutron Capture Therapy (BNCT) Monte. Carlo treatment simulations were investigated. Presently, the Idaho National Engineering and Environmental Laboratory (INEEL) and Montana State University (MSU) treatment planning software uses Non-Uniform Rational B-Spline (NURBS) surfaces to allow construction of a patient model from medical image scans. Radiation doses are then calculated using a Monte Carlo transport method that employs a ray-tracing algorithm to find intersections in the NURBS geometry. To obtain statistically significant results, a Monte Carlo treatment simulation can easily require several hours. A profile of the simulation has revealed that between 65 and 90 percent of the time is spent in computing particle intersections with NURBS surfaces. Performance tuning in this area could result in substantial time savings. The new algorithms track rays through uniform volume elements (univels) rather than finding intersections with NURBS surfaces. Line rasterization methods have been adapted to allow fast tracking through the univels. Simulation times have been measured on Pentium Pro systems with greater than 64 MB of RAM. Resultant accuracy was measured by comparing the results to that of a problem with an accepted answer. Results indicate that the new univel-based tracking methods are approximately 50 times faster than the NURBS-based methods. By incorporating these new tracking methods into the Monte Carlo simulations, speedups by factors of 6 to 10 are common. Accuracy was maintained by scaling the univels to match the precision of the original image data. The results of the new simulations have been shown to be well within the accepted experimental error. The new univel-based particle tracking methods outperform the NURBS-based methods. The univel geometry is a simple representation that can be easily generated. The freedom of independently labeling individual univels as a particular material allows creation of arbitrarily shaped bodies not easily represented by a NURBS surface. Because of the speed and success of the newer tracking methods, the future BNCT planning process will be entirely univel-based. Many of the developed tools, methods, and representations could also be applied to non-medical problems and other transport simulations.  RAPID GEOMETRY INTERROGATION FOR A UNIFORM VOLUME ELEMENT- BASED MONTE CARLO PARTICLE TRANSPORT SIMULATION by Michael William Frandsen A thesis submitted in partial fulfillment of the requirements for the degree of Master of Science in Computer Science MONTANA STATE UNIVERSITY—BOZEMAN Bozeman, Montana April 1998 K31V F^6 ii APPROVAL of a thesis submitted by Michael William Frandsen This thesis has been read by each member of the thesis committee and has been found to be satisfactory regarding content, English usage, format, citations, bibliographic style, and consistency, and is ready for submission to the College of Graduate Studies. Dr. Gary H. Harkin Approved for the Department of Computer Science Dr. J. Denbigh Starkey Approved for the College of Graduate Studies Dr. Joseph J. Fedock STATEMENT OF PERMISSION TO USE In presenting this thesis in partial fulfillment of the requirements for a master's degree at Montana State University—Bozeman, I agree that the Library shall make it available to borrowers under rules of the Library. IfI have indicated my intention to copyright this thesis by including a copyright notice page, copying is allowable only for scholarly purposes, consistent with "fair use" as prescribed in the U.S. Copyright Law. Requests for permission for extended quotation from or reproduction of this thesis in whole or in part may be granted only by the copyright holder. Signature ^ DateIZ iA m ________ ACKNOWLEDGEMENTS This work was sponsored by the U.S. D.O.E., Office of Energy Research, through INEEL under Contract DE-AC07-94ID13223. Thanks are given to Floyd Wheeler of the INEEL for providing the Monte Carlo transport simulation results and to Dan Wessol, also of the INEEL, for introducing me to the problem at hand and starting me down the path of uniform volume element tracking. Thanks are also given to my wife for her constant support and assistance in editing this thesis - not to mention the many hours she let me spend in front of her computer instead of with her. TABLE OF CONTENTS Page BACKGROUND INFORMATION....................................................................................I Boron Neutron Capture....................................................................................................I Boron Neutron Capture Therapy.....................................................................................I Boron Neutron Capture Therapy for Glioblastoma Multifbrme.................................... 2 Clinical Trials....................................................... :............................................... .........3 BNCT Planning Process.......................................................................... i..................... 3 BNCT Treatment Plan Validation.......................................................................... 4 INTRODUCTION............................................................ 5 Problem Statement.......................................................................................................... 5 ANALYSIS OF CURRENT TRACKING TECHNIQUE................................................. 7 Non-Uniform Rational B-Spline Surface Descriptions.................................................. 7 Positive Aspects of NURBS.................................................................... 7 Problems with NURBS................................................................................................... 8 UNIFORM VOLUME ELEMENTS PROPOSED...........................................................10 Description.....................................................................................................................10 Advantages....................................................................................................................11 Disadvantages................................................................................................................13 TWO-DIMENSIONAL TRACKING TECHNIQUES........................................... 14 Common Integer-Based Representative Line Algorithm.............................................. 15 Adaptation to Investigate All Pixels Intersected by Ideal Line..... ............................... 17 Relation to Ray Tracking.............................................................................................. 18 Discussion of Extension to Non-Grid Points.................................................................19 Discussion of Two-Dimensional Tracking Error.......................................................... 20 EXTENSION TO TRACKING IN THREE DIMENSIONS........................................... 22 Sample Problem........................................................... 22 , Representative Elements Intersected by Ideal Line Examined..................................... 23 Representative Elements Intersected for Sample Problem........................................... 25 All Elements Intersected by Ideal Line Examined........................................................ 27 All Elements Intersected for Sample Problem........ ,.................................................... 31 RESULTS...................................................................................................................... ...33 How the Algorithms Fared........................................................................................ i-..33 ' Error Analysis of Only Examining Representative Univels along the Ideal Line........ 39 Analysis of Greater Distance Covered in Examining Only Representative Univels Along the Ideal Line.............. 40 Monte Carlo Transport Simulation Results.................................................................. 42 CONCLUSIONS.............................................................................................................. 46 Future Directions.......................................................................................................... 47 REFERENCES................................................................................................................. 48 APPENDIX....................................... 49 Three-Dimensional Algorithm to Examine Representative Univels along Line.......... 49 Three-Dimensional Algorithm to Examine All Univels Intersected by Line................51 Timing Comparisons of Various Algorithms through Same Number of Elements...... 54 vi vii LIST OF TABLES Page Table I : Univel Locations and Sampling Points for Sample Problem.............................25 Table 2: Univel Locations, Entrant Face Intersections, and Distance Inside each Univel for Sample Problem..................................................................................................31 Table 3: Best to Worst Performance for 1000 Univels.................................................... 38 Table 4: Best to Worst Performance for Initialization..................................................... 38 Table 5: Best to Worst Performance for 10 Univels........................................................ 38 Table 6: Best to Worst Performance for 20 Univels........................................................ 38 Table 7: Processor Time Requirements for Patient Irradiation Simulation......................44 Viii LIST OF FIGURES Page Figure I : Uniform Volume Element Structure and Storage.............................................10 Figure 2: From Medical Images to Univels............................... ...........*.......................... 12 Figure 3: Sides "Owned" by a Pixel and a Univel.............................................................15 Figure 4: Two Possibilities for Drawing Pixels along a Line............................................15 Figure 5: Integer-Based Algorithm to Plot Representative Line between Grid Points..... 16 Figure 6: Integer-Based Algorithm to Plot All Pixels along Line between Grid Points....18 Figure 7: Maximum Skipped Segment Length within a Pixel.......................................... 21 Figure 8: Sample Three-Dimensional Line....................................................................... 22 Figure 9: Calculation of an Initial Error Term for the Centered Y-Value........................ 24 Figure 10: Graphical View of Intersected Representative Univels for Sample Problem ..26 Figure 11: Distances to Each Candidate Face.................................................................. 28 Figure 12: Examination of Error Terms to find which Coordinate to Increment............. 29 Figure 13 : Graphical View of All Intersected Univels for Sample Problem.................... 32 Figure 14: Long-Term Performance of Various Tracking Algorithms............................ 35 Figure 15: Typical Performance of Candidate Tracking Algorithms............................... 36 Figure 16: Wireframe Representation of Surfaces for Head, Brain and Target Regions ..43 Figure 17: Two-Dimensional Slice Showing Relationship of BMRR Beam Delimiter for the rtt_MC Model of the Head with an 8-cm Depth Target.....................................43 ABSTRACT New particle tracking methods to speed Boron Neutron Capture Therapy (BNCT) Monte. Carlo treatment simulations were investigated. Presently, the Idaho National Engineering and Environmental Laboratory (INEEL) and Montana State University (MSU) treatment planning software uses Non-Uniform Rational B-Spline (NURBS) surfaces to allow construction of a patient model from medical image scans. Radiation doses are then calculated using a Monte Carlo transport method that employs a ray­ tracing algorithm to find intersections in the NURBS geometry. To obtain statistically significant results, a Monte Carlo treatment simulation can easily require several hours. A profile of the simulation has revealed that between 65 and 90 percent of the time is spent in computing particle intersections with NURBS surfaces. Performance tuning in this area could result in substantial time savings. The new algorithms track rays through uniform volume elements (univels) rather than finding intersections with NURBS surfaces. Line rasterization methods have been adapted to allow fast tracking through the univels. Simulation times have been measured on Pentium Pro systems with greater than 64 MB of RAM. Resultant accuracy was measured by comparing the results to that of a problem with an accepted answer. Results indicate that the new univel-based tracking methods are approximately 50 times faster than the NURBS-based methods. By incorporating these new tracking methods into the Monte Carlo simulations, speedups by factors of 6 to 10 are common. Accuracy was maintained by scaling the univels to match the precision of the original image data. The results of the new simulations have been shown to be well within the accepted experimental error. The new univel-based particle tracking methods outperform the NURBS-based methods. The univel geometry is a simple representation that can be easily generated. The freedom of independently labeling individual univels as a particular material allows creation of arbitrarily shaped bodies not easily represented by a NURBS surface. Because of the speed and success of the newer tracking methods, the future BNCT planning process will be entirely univel-based. Many of the developed tools, methods, and representations could also be applied to non-medical problems and other transport simulations. I BACKGROUND INFORMATION Boron Neutron Capture The principle of Boron Neutron Capture (BNC) is that a neutron striking a 10B nucleus can produce 11B which fissions into Lithium (7Li) and Helium (a-particle) nuclei. The probability of neutron capture by 10B is orders of magnitude greater than for constituent nuclides of biological matter and this probability is inversely proportional to the neutron kinetic energy of the neutron. The Lithium and Helium nuclei are directed opposite each other with significant energy. Boron Neutron Capture Therapy Boron Neutron Capture Therapy (BNCT) is a method that uses this reaction as a treatment modality. 10B is introduced to the patient by injection of an appropriate amount of one of a variety of Boron compounds that are not toxic for prescribed doses. After some time has elapsed, the Boron compound has been delivered throughout the body through the bloodstream. The compound is engineered to target the tumor cells. Also, because cancer cells have an abnormally high growth rate, they take in more of the Boron compound than nearby healthy tissue. To produce the reactions, a neutron beam is directed at the cancer area. The treatment consists of turning the beam on for an amount of time, typically one-half to one hour, determined by the beam’s neutron intensity. The Boron Neutron Capture reaction inactivates the host cell and, with much lower probability, some nearby cells when the Boron fissions. The beam path and concentration 2 of the Boron compound in the cancer cells focuses on destroying cancerous tumor cells and sparing healthy tissue to the greatest extent possible. Boron Neutron Capture Therapy for Glioblastoma Multiforme Montana State University (MSU) and the Idaho National Engineering and Environmental Laboratory (INEEL) have been involved with creating the user interface to the planning software for treating patients with a specific brain tumor called Glioblastoma Multiforme since 1991. This cancer affects 7000 people in the U.S. each year and is generally considered fatal within six months if untreated or perhaps one year even with treatment. The reason Boron Neutron Capture Therapy is a treatment candidate for Glioblastoma Multiforme is multi-faceted. Conventional methods of therapy are not optimal. Removal of the tumor through surgery alone often results in tumor regrowth due to the difficulty of removing the small tumor tendrils in nearby tissue without also removing much of the surrounding healthy brain mass. In Boron Neutron Capture Therapy, some researchers consider the neutron beam as effective on the tumor tendrils as it is on the tumor bulk although the bulk is often removed through surgery prior to Boron Neutron Capture Therapy. The Boron compound concentrates in the tumor and therefore the cell kill is concentrated only in areas near the tumor. Only cells within the range of the reaction products, approximately one cell diameter, are destroyed. 3 Clinical Trials The early-stage clinical trials of Boron Neutron Capture Therapy for treatment of Glioblastoma Multiforme are promising. At this time, clinical trials have been conducted on only 36 patients at Brookhaven National Laboratory (BNL) during the initial U. S. clinical trials. In the future, there may be a much greater demand for this type of treatment. Additionally, the method is also being explored for treating or assisting in the treatment of other types of cancer such as lung cancer. This potential rise in demand for Boron Neutron Capture Therapy will require planning software optimized for patient throughput. Optimizations to the current software are necessary to meet these demands. BNCT Planning Process Although Boron Neutron Capture Therapy itself can be conducted in a short session, the planning for that treatment has many stages. The patient must undergo a medical scan such as Nuclear Magnetic Resonance Imaging (NMRI) of Computed Tomography (CT). Next, physicians determine the locations of important structures such as the scalp, skull, brain, target area, and tumor by analyzing the digitized medical scans. Additionally, critical areas are marked for avoidance. The physician often does the marking by outlining structures of interest on several medical images. The stacked slices are then representative of the patient's head and internal regions of interest. Three- dimensional visualization software can be used to observe and verify the location of the defined structures. The three-dimensional views are also useful for the physician to 4 determine an appropriate beam intensity and beam direction for the therapy. BNCT Treatment Plan Validation Before the treatment proceeds, simulation software must determine the optimal irradiation geometry and irradiation time. A Monte Carlo simulation developed at the Idaho National Engineering and Environmental Laboratory is used to track millions of particle paths and calculate the radiation doses associated with the treatment. Radiation isodoses are then superimposed over the original medical images to aid physicians in development of a treatment plan. Although the planning software interface developed by Montana State University and the Idaho National Engineering and Environmental Laboratory is widely used throughout the world for these types of treatments, enhancements could be made to minimize the time required for the planning process. 5 INTRODUCTION The method presented here is used primarily to model neutron/photon transport for BNCT although it could be applied to model most transport processes where the geometry is first decomposed into univels. The INEEL and MSU BNCT collaboration has developed a graphical environment for reconstructing univel and NURBS-based solid models derived from medical images. These solid models form the geometric basis of the tumor-target areas and surrounding healthy tissue for the BNCT [1] treatment planning radiation transport model [2,3]. The system manipulates and samples medical images so that accurate geometric models of imaged anatomical features can be created and viewed. Problem Statement Computer simulations that model a proposed treatment plan to determine its efficacy are part of the planning stage of BNCT for glioblastoma treatment. It has been found that these simulations can take substantial amounts of time. One major component of these simulations is the need to track simulated rays through the materials of interest and report intersection distances to these regions. A typical simulation on a high-speed computer may require three hours to track an adequate number of particles. It has been determined that 65% to 90% of this time is spent in tracking particles from one region of interest in the NURBS geometry to the next and returning a distance to this intersection. If this bottleneck can be reduced, the overall simulation time can be reduced dramatically. This thesis is an investigation and construction of techniques to decrease the computational time required by Boron Neutron Capture Therapy Monte Carlo transport simulations by optimizing the particle tracking routines. 6 7 ANALYSIS OF CURRENT TRACKING TECHNIQUE Non-Uniform Rational B-Spline Surface Descriptions The current tracking routine is based on interrogating Non-Uniform Rational B- Spline (NURBS) surfaces. Using NURBS is one method of fitting a smooth curve to a set of control points. Many NURBS curves at different heights can be joined to produce a NURBS surface. Only a relatively small set of control points is required to describe a NURBS surface. For BNCT, the NURBS control points for a given structure are generally defined on 5 to 40 parallel medical images depending on the size of the structure. Each object in each slice may be defined with 10 to 20 control points, yielding approximately 50 to 800 control points to describe a given structure. The number of control points will depend on both the size of the structure and how accurately the structure must be represented. Given these control points and a degree for the NURBS surface, a resultant outer boundary of the structure is defined. Generally, a quadratic fit is used for the control points. A linear fit would produce a flat-faced polygon whereas a cubic or higher fit does not follow the control points as closely. Therefore, a quadratic fit seems most representative of the underlying surfaces. Positive Aspects of NURBS For many reasons, NURBS are an appropriate representation for the surfaces of the regions of interest. Even with a small number of control points, the NURBS surface 8 can fairly accurately depict most anatomical structures of interest. The NURBS exhibit very good local control and can be easily adjusted with the control points. The resultant structures lend themselves well to smooth three-dimensional surface renderings for visualization. The compact representation consisting of a small number of ordered control points and degree for the NURBS surface requires little memory to describe the geometry. Problems with NURBS Until the last few years, computer memory has been at a premium. Another reason for the original selection of NURBS surfaces is due to the greater memory requirements of many other representations. Today, this is not nearly as much of a problem. On the contrary, users today suffer because of the trade-off that storage is minimized at the expense of requiring more computation time when the surface is computed or the geometry is interrogated. The control points are a minimum description for the NURBS mathematical structure. A lot of computations are required to determine if a given point is inside, outside, or on a given NURBS surface. Current methods used for determining the first point of intersection of a ray with a NURBS surface use bounding boxes of decreasing size to either hone in on the point of intersection or determine that no intersection occurs. This method of computing intersections is valid and has been used for the treatment simulations but is time-consuming and not without some other flaws. In addition to the high computation requirements of computing intersections with NURBS surfaces there are problems of lost particles and overlapping structures. Due to 9 precision problems, lost particles are common in NURBS geometries when surfaces are near each other. Somewhat related to this is the problem of overlapping surfaces. A user could inadvertently create control points that cause close NURBS surfaces to overlap. When overlap occurs, the volume within the overlap is indeterminate as to what structure it belongs to. Another problem is that certain types of structures are difficult to represent with NURBS. Structures that branch and structures with multiple connected components produce problems. Typically, these surfaces are either approximated by something that the NURBS can describe or require multiple NURBS for their definition. If multiple NURBS are used, gaps or overlaps could be present where the underlying structures should be connected. Given these problems with both the interrogation time requirements and potential ambiguities in the description, it is worthwhile to investigate other representations for the regions of interest and associated tracking methods. 10 UNIFORM VOLUME ELEMENTS PROPOSED Description A uniform volume element (univel) will describe the volume enclosed in a right parallelepiped (or rectangular box). The term voxel is not used because just as pixels are generally considered squares, voxels are generally considered cubes. Here, each uniform volume element has the same length, width, and height but the shape is not necessarily a cube. Many univels can be placed adjacent to each other to describe a volume. Even though each univel is box-like, the elements can describe a structure of arbitrary precision determined by the univel size. Greater accuracy will require more storage. Figure I shows what a typical univel might look like along with its normalized representation for array storage. z A ẐTl c mm , b mm I/ a mm (a) Univel > / / Z I (b) Array storage for univel Figure I: Uniform Volume Element Structure and Storage 11 Advantages Univels are proposed as an alternative to the NURBS representation because of the following advantages. The univel representation of a volume can easily be thought of as a three-dimensional array of volume elements. In this array representation, it is computationally simple to find what structure is present at a given point in space. Intersections of a ray with an edge of a structure would require tracking the ray through the three-dimensional array until a volume element of a new material type is encountered. The array mathematics is considerably simpler and also quicker than the NURBS mathematics required for tracking an equivalent ray. The univel representation can also be made to have a simple correspondence with the original medical image slices. If the univels are to be the same resolution as the medical image slices, each pixel on a medical image slice will correspond directly to one univel. The univels can then be stored and retrieved much the same as the original medical image slices. Using univels of finer resolution is certainly possible but the structure would be determined more accurately than the original data. The amount of data available does not lend itself well to support this gain in the resolution of the resultant defined structures. Using univels of a lesser resolution would reduce geometric accuracy. Another advantage is that each univel can be labeled independently as any particular material. There cannot be a problem of overlapping regions and there is no problem of adjusting control points for a surface boundary to accurately fit an underlying structure. Rather than using control points, each individual element is labeled with an 12 appropriate material. Finally, the memory demands of the univels are approximately the same as the original medical images and by no means a problem on today's computer systems. Figure 2 illustrates how the univels relate to the original medical images. In this case, Figure 2a shows three slices each of 4x4 resolution that comprise the stacked image set. For each pixel, a corresponding univel will be created. Figure 2b shows what the stacked univels might look like. The dashed lines indicate the positions of the original medical image slices that are usually located at the transverse midpoint of the univels. Figure 2c illustrates a graphical representation of the three-dimensional array that the univel information will be stored in. (a) Original medical (b) Stacked univels image slices / / / / / / / / / / / / / / / Z / / / / / Z / / / Z / Z Z Iz Z Iz (c) Array storage for univels Figure 2: From Medical Images to Univels 13 Disadvantages Although univels are useful for the purposes described herein, they do have their disadvantages. The storage required for univels is greater than that of many other volume (or surface) representations. The box-like shape does not preserve the surface area of the represented volume. Additionally, curved surfaces are stair-stepped to a degree depending on the size of the univels. Because of the stair-steps, computation of normal vectors is not very informative. However, these are not concerns for the Monte Carlo simulations. 14 TWO-DIMENSIONAL TRACKING TECHNIQUES It is simpler to first analyze the two-dimensional problem before discussing tracking in three dimensions. The two-dimensional problem is similar to drawing a line of pixels on a computer graphics display. Because this is such a common task, many algorithms were developed for efficiently drawing lines. Some references consider a pixel as a grid point occurring at the intersection of two gridlines. Other references consider a pixel to cover the area between two adjacent horizontal and vertical gridlines. Here, the latter view will be taken. Assuming the grid is composed of horizontal and vertical lines for each integral value of x and y, the center of any pixel lies at (1+1/2, J+l/2) where I and J are arbitraiy integers. However, the entire pixel can simply be referenced by (I, J). To settle the ambiguity of which edges belong to a pixel and which to its adjacent neighbors, some heuristic must be used. As recommended by the Graphics Gems series [4], the floor function will be used to determine which pixel (or later which volume element) contains a given floating point ordered pair (or ordered triple). In other words, the pixel denoted by (I, J) will span the area enclosed by: I<=x0 Eliminate y errxz<0 Eliminate z erryz<0 Eliminate z errxz>0 Eliminate y errxy<0 Eliminate y errxy>0 Eliminate x Figure 12: Examination of Error Terms to find which Coordinate to Increment 30 In the algorithm, one error is arbitrarily checked first and the sign determines which error term is to be checked next. The coordinate to update is determined after just two comparisons. As this coordinate is stepped by one unit, it has an effect on two of the three error terms. The error terms are updated to indicate a step was taken in the calculated direction. As before, not much attention is paid as to what should be done if any error is precisely zero, indicating that more than one component should be increased. Here, speed is gained by using fewer comparisons for the trade-off that a rare univel will be investigated that only touched the line at a single point. However, it is not difficult to explicitly test the case of the errors being zero and incrementing more than one component as necessary. One aspect of this algorithm can be exploited that is not of use in the case of traversal or simply marking examined univels. Since the traversal is ultimately used to determine an intersection point, knowing which face of the univel was first intersected can save time. Here, that can be easily managed. The algorithm that visits representative univels will require looking at each of three candidate faces and will therefore have some additional overhead when determining the final intersection. Although this algorithm tracks through more univels for a given segment length, the final calculation of an intersection point can be computed more quickly. Code for the developed algorithm is given in the appendix. 31 All Elements Intersected for Sample Problem For the sample problem depicted in Figure 8, Table 2 lists the array positions that will be examined and the actual coordinates of the face intersections in the case of examining all univels intersected by the ideal line. Additionally, the length of the line segment inside each univel is given. Table 2: Univel Locations, Entrant Face Intersections, and Distance Inside each Univel for Sample Problem Arrayx Arrayy Arrayz X Y Z Distance 0 0 2 — — — 0.199 0 I 2 0.056 1.000 2.594 0.863 0 I 3 0.296 1.722 3.000 0.332 0 2 3 0.389 2.000 3.156 1.195 0 3 3 0.722 3.000 3.719 0.597 0 3 4 0.889 3.500 . 4.000 0.398 I 3 4 1.000 3.833 4.188 0.199 I 4 4 1.056 4.000 4.281 1.195 I 5 4 1.389 5.000 4.844 0.332 I 5 5 1.481 5.278 5.000 0.863 I 6 5 1.722 6.000 5.406 0.996 2 6 5 2.000 6.833 5.875 0.199 32 All univels examined for the sample problem are depicted in Figure 13. £__y Figure 13: Graphical View of All Intersected Univels for Sample Problem 33 RESULTS The results will be divided into two main parts. In the first part, various tracking algorithms will be compared. The algorithm that investigates representative univels along the ideal line was determined to be most appropriate for the Monte Carlo simulations. In the second part, the results of using this new technique in the Monte Carlo transport simulations will be discussed. How the Algorithms Fared Time trials were conducted for the described algorithms as well as others. The algorithms are classified according to four categories. Two-dimensional algorithms will be designated as 2D and three-dimensional algorithms as 3D. Although grid point-based algorithms are not accurate enough for the tracking, they were used for comparison. T will indicate that the algorithm started at an integer grid point with integer increments and ‘F ’ will be used for an arbitrary starting point in an arbitrary direction. A second T or ‘F ’ will indicate if the algorithm uses integer arithmetic or floating point arithmetic in its inner loop. Finally, ‘A’ will indicate that the algorithm examines all elements along the path and R will indicate only representative elements were examined. Hence, 3DFIR will indicate the three-dimensional algorithm that starts at an arbitrary position, is directed in an arbitrary direction, uses integer arithmetic for its main stepping loop, and only tracks through representative elements. The integer starting point algorithms will be used as a basis of comparison for the 34 floating start point algorithms. If times are close then it can be concluded that the generalized algorithms cannot perform much better. This assumes that the entirely integer-based algorithms are very near a lower bound for this type of tracking environment. Tracking through NURBS will not be analyzed in this section since it will be evident in the next section that the univel tracking methods are much faster. However, the question for the time being will be: How much of a performance trade-off do the algorithms take for the necessary added starting point precision and tracking precision? The tests are a somewhat simplified version of shooting many random rays. To avoid undue complication, the algorithms used assume that x is the dominant direction of the ray and that each component of the direction vector is positive. It is assumed that through symmetry the other cases need not be considered. Another simplification is that the algorithms do not calculate a precise intersection with the final univel on the path. It is the tracking that is the greatest concern here and so it is assumed that a small amount of constant time can be added to each algorithm to account for the calculation of the final intersection point. Of course, this could vary slightly between the algorithms. The test was conducted by generating a random data set of qualifying starting points and direction vectors. For the pure integer versions, the values of the starting points were truncated to integers and the value of each component of the direction vector was multiplied by a large constant and each component then truncated. The number of pixels (in the two-dimensional case) and univels (in the three-dimensional case) was varied over a large range of values and corresponding times required by each respective algorithm were recorded. The tests were run on a Linux-based Pentium Pro 150 with 96 35 MB of RAM. A plot of all times showing long-term performance is given in Figure 14. Tracking Times through Large Number of Elements, PPRO 150 50.00 40.00 30.00 20.00 10.00 Elements Examined ---B--- 2D MR - - -B - - - 2D NA --a-- 3D MR 3D FIR 3D FFR 3D MA —•— 3D FIA 3D FFA Figure 14: Long-Term Performance of Various Tracking Algorithms 36 Additionally, a plot comparing the candidate algorithms (3DF—) for a smaller number of univels is given in Figure 15. Tracking Times through Small Number of Univels, PPRO 150 Univels Examined —G— 3D FIR 3D FFR -■-3D FIA — 3D FFA Figure 15: Typical Performance of Candidate Tracking Algorithms I 37 There is a comparison problem that examining an equal number of univels does not mean equal distances were covered. When only representative univels are examined, a step from one element to the next is typically larger than if all stepped through univels are face-adjacent to each other. A correction factor for this will be examined later. The tests as given are intended to show how the algorithms compare for the same number of elements or, equivalently, for the same number of iterations. As expected, the long-term performance of the four algorithms that relied on floating point calculations in their main loop performed worst. The algorithms that used integer arithmetic in their main loops but added the precision of using floating point starting points and floating point direction vectors performed, on average, approximately the same as their counterparts that used integers throughout. The y-intercept of each line indicates the initialization time required for each algorithm. This is of importance when only a small number of univels are to be examined for a given intersection test. Fifteen univels is a good estimate of an average number of univels that must be examined for a given intersection test in the current Monte Carlo simulations. Tables 3,4, 5, and 6 summarize the performance of four candidate algorithms in the cases of long-term performance, initialization performance, performance for 10 univels, and performance for 20 univels, respectively. 38 Table 3: Best to Worst Performance for 1000 Univels 3DFIR 3DFIA 3DFFA 3DFFR Relative Time 1.00 1.32 3.93 7.18 Table 4: Best to Worst Performance for Initialization 3DFFA 3DFFR 3DFIA 3DFIR Relative Time 1.00 1.42 1.65 2.12 Table 5: Best to Worst Performance for 10 Univels 3DFIA 3DFIR 3DFFA 3DFFR Relative Time 1.00 1.14 1.25 2.10 Table 6: Best to Worst Performance for 20 Univels 3DFIA 3DFIR 3DFFA 3DFFR Relative Time 1.00 1.04 1.62 2.80 39 Somewhere between 3 to 7 univels is where both 3DFI algorithms start having better performance than the two 3DFF algorithms. Since it is assumed that most intersections will require investigating more than 3 to 7 univels and since the longer intersections will have more of an effect on the overall time, it is evident that the algorithms that use integer arithmetic in their inner loops should be used in most cases. For the general problem, 3DFIA and 3DFIR are both candidates. Of course, 3DFIR can only be used when it is allowable to skip some univels and miss some potential intersections. Rates given so far measure univels examined. Another measure is total tracking distance. 3DFIA does not travel as far as 3DFIR when tracking through the same number of univels. A multiplication factor will be derived to compare these algorithms on equal terms of distance covered. More data for the timing tests are given in the appendix. Error Analysis of Only Examining Representative Univels along the Ideal Line For the Monte Carlo problem at hand, intersections missed by 3DFIR have been found to not be a problem. On the contrary, they have more good side effects than bad. In the case of tracking through regions of the head, locality of regions is assumed. That is, regions are generally found in localized globs rather than random spread-out points. In this case, when a new material is entered, it is typical that the continuation of this line will continue to pass through at least a few univels of the same material. For this reason, 3DFIR does not miss most intersections although it may appear that a missed univel could result in a slightly greater distance-to-intersection reported. However, three 40 candidate planes are examined when the point of intersection is being determined. The candidate planes do not necessarily correspond to the calculated univel. They are simply the past three x-y, x-z, and y-z planes, or faces, intersected by the line and so correspond to the current univel and two previous univels (not necessarily the past two intersected univels). This examination of previous univels makes it likely that the resultant intersection point will be exact. When an intersection is missed entirely, it is due to the ray entering the new material only briefly through a small portion of a univel and then immediately returning to the original material. For a univel of dimensions a mm x b mm x c mm, it is relatively straightforward to calculate that all univels intersected by a line segment of at least sqrt(a2+b2+c2)/2 mm are examined. This is equal to half the length of the main diagonal of the univel. See Figure 7 for the two-dimensional analog. In the case of the missed intersections, the amount of material through which the ray passes is so small that it is not necessary to proceed with the transport calculations for the intersection. On the contrary, these extra calculations do little to affect the overall problem but will increase the time required for the simulation. This, of course, relies on a small univel size compared to the mean-free path of the particles being examined. This is satisfied in the current simulations and should hold true in the future assuming imaging devices remain unchanged or become more resolute. Analysis of Greater Distance Covered in Examining Only Representative Univels Along the Ideal Line Earlier comparisons gave timings based on different algorithms tracking through 41 the same number of elements. However, tracking through only representative elements rather than all elements covers more distance for the same number of univels examined. In general, however, a simple multiplication factor can be derived to estimate the proportionality constant that exists between the two distances covered for the same number of univels. It is easy to analyze the extreme cases when the ray moves parallel to an axis or when the ray moves equally in all directions. If the ray moves parallel to an axis, both algorithms track through all univels and the proportionality constant is therefore 1.0. In the case of all directions varying equally, it takes three steps for the algorithm tracking through all univels to increase each of x, y, and z The algorithm investigating representative elements will increase whichever component it arbitrarily decided was changing most rapidly and at the same time see that it must also correspondingly increase each of the other two components. What was achieved in one step by the representative tracking algorithm took three steps in the other algorithm. Hence, the representative tracking algorithm covers a distance 3.0 times as great as the competing algorithm in this worst case while stepping through the same number of univels. Therefore, the average proportionality constant must lie between these two extremes of 1.0 and 3.0. To think of this another way, this is similar to considering two different metrics for distance from a source univel to a destination univel. The algorithm tracking through all elements takes the unit distance between univels in each of x, y, and z and sums them to calculate the number of univels examined. The representative tracking algorithm examines only the maximum of the three components. 42 For the continuous case, assume a sphere of radius r. The continuous and discrete cases should be similar Assume that a ray starts at the center of the sphere and moves in an arbitrary direction to its surface. The ratio between the average number of univels tracked through in 3DFIA and 3DFIR is: averageflxl+lvl+lzn over the surface of the sphere average(max(|x|, |y|, |z|)) over the surface of the sphere In both cases, x2+y2+z2=r2. Empirical trials show that this ratio is approximately 1.95. In other words, the representative element tracking ray will travel about twice as far, on average, while investigating the same number of univels as the algorithm that tracks through all univels. Because of this substantial difference, it is optimal to use the representative tracking algorithm whenever it is estimated that more than 4 univels will be tracked through in computing an intersection. (Tracking through 4 univels with 3DFIR takes approximately the same time as tracking through 8 univels with 3DFIA based on the appendix data.) However, the difference is not substantial when tracking through a small number of univels. If tracking through 10 univels in the 3DFIR algorithm requires I time unit, equivalent tracking through 20 univels in the 3DFIA algorithm requires approximately 1.17 time units. Monte Carlo Transport Simulation Results The univel geometry and representative univel tracking technique, 3DFIR, was tested using the rttJMC Monte Carlo module [8] of the BNCT treatment planning software. Testing was performed using Linux-based Pentium Pro systems. A model, 43 previously constructed from medical images of a male volunteer, was used to compare processor times and accuracy in the simulated irradiation. Figure 16 shows a wire frame representation of the NURBS surfaces originally reconstructed from the images. For this simple model only the outer surfaces of the head, the brain, and a simulated target volume were represented. Figure 17 shows a two-dimensional raster image of the model including the neutron beam delimiter in use for the human clinical trials at Brookhaven National Laboratory Medical Research Reactor (BMRR) [9]. Figure 16: Wireframe Representation of Surfaces for Head, Brain and Target Regions Figure 17: Two-Dimensional Slice Showing Relationship of BMRR Beam Delimiter for the rtt MC Model of the Head with an 8-cm Depth Target 44 Three different geometric representations were calculated. For the current patient calculations, the NURBS representation is used for the patient geometry and analytical surfaces. Combinatorial Geometry (CG) [10] is used to represent the neutron beam regions. These results were compared to a case where just the patient geometry used the univel representation and to a case where all regions, including the beam geometry, were represented as univels. For each case, 2x106 neutron history simulations and IxlO6 gamma ray simulations were performed. Timing comparisons are shown in Table 7. Table 7: Processor Time Requirements for Patient Irradiation Simulation Geometry Patient Representation Beam Time (Minutes) Time Relative to NURBS NURBS CG 332 1/1 Univels CG 41 1/8 Univels Univels 22 1/15 In some cases, due to round off, the Monte Carlo software detects an impossible particle location situation and must discard the particle since it is not possible to know what to do next. This is referred to as a Tost particle’ and Tost particles’ are common in NURBS geometries when surfaces are near each other. For this simulation, there were four Tost particles’ for the NURBS case and none for the univel cases. Often, there are more Tost particles’ in NURBS geometries. Lost particles have yet to be encountered using the univel tracking routines. This indicates, with all things otherwise being equal, 45 that the univel method will be more precise. For these three cases, a detailed check showed all results to be within the statistical errors of the Monte Carlo process. A profile of the time requirements showed that the relative ray trace time for the NURBS geometry was 64%. This was reduced to 6% for the case with univels. In another test case, a timing comparison was made using a different number of univels to determine the sensitivity of time requirements to the number of univels, or precision of the geometric representation. For this case, the geometry was an H2O cylinder with a simulated irradiation in the Brookhaven neutron beam. For this simple geometry problem, the CG model required about 20% less time than the univel model using 13,520 univels. When the spacing along the cylinder axis was reduced to construct a univel model having 3.28x106 univels (a factor of 243 increase in the number of univels), the processing time was increased by only 74%. For more complicated models where several CG objects are required to represent the geometry, univel tracking will be more efficient than CG tracking. 46 CONCLUSIONS Univel-based modeling relying on uniformity and integer-based arithmetic works well and tracking techniques based on this are substantially faster than the methods based on more compact NURBS representations that rely heavily on floating point arithmetic. Because of the significant gain in speed with no significant increase in error, the representative univel tracking method is a good method for Monte Carlo transport simulations. Skipping over some univels has little effect on accuracy. Those times that intersections were missed, it was probably better to skip them as they proved to be insignificant and would only increase computation time with little or no gain in accuracy. A typical I mm resolution is already quite fine and so skipping such a tiny volume element is not substantial, especially considering that the original segmentation of the images cannot be expected to be accurate to the original pixel level. Tables I and 2 show that the lengths of line segments through the skipped univels are short compared to the lengths through the examined univels. In the case that greater accuracy is obtained in defining the regions, the smaller univels will allow increased accuracy at the expense of requiring more space and slightly more time. With the new methods, such a small fraction of time is spent in ray tracking that further optimizations would do little to affect the simulation time. This would only be worthwhile if the number of univels used to represent the geometry grows by a substantial factor. If 20 univels were tracked through on average at a 256x256x40 resolution, 80 univels would be tracked through on average at a 1024x1024x160 resolution. Even with this tremendous jump, the amount of time spent tracking would barely double. The overall simulation would hardly be affected 47 since tracking has now been reduced to approximately 6% of the total simulation time. Future Directions Using a derivative of the 3DFIA or 3DFFR algorithms may optimize a remaining portion of the Monte Carlo simulations. The greatest percentage of simulation time is now spent in the tally portion of the code. One part of this Stage requires voxels to be stepped through, one at a time, calculating the distance through each voxel. Both 3DFIA and 3DFFR algorithms are ideal for this type of tracking since they are fast, track through all voxels, and can quickly compute distances through voxels by keeping track of entered and exited faces. Use of these algorithms for the tally portions of the simulations is being investigated and looks promising, In the current simulations, only a small number of voxels are tracked through so the speedup may be only a small factor. However, large speedups may be anticipated in future simulations that may require deeper penetrating particles. 48 REFERENCES 1. Wessol, D.E.; Wheeler, FJ.; Babcock, R.S; Starkey, J.D., “Medical Image Reconstruction and Geometry Generation Software”, 5th International Symposium on Neutron Capture Therapy, Columbus, Ohio, September 13-17, 1992. 2. Wessol, D.E.; Babcock, R.S.; Wheeler, F.J.; Harkin, GJ.; Voss, L.L.; Frandsen, M.W., “BNCT_Rtpe: BNCT Radiation Treatment Planning Environment User’s Manual”, Version 2.2, February 21,1997, http://id.inel.gov/cart/rtpe- manual/secOO.html. 3. Capala, J.; Chadha, M.; Wheeler, FJ.; Wessol, D.E.; Atkinson, CA., “Radiotherapy Planning for the Brookhaven Trials of BNCT for Human Gliomas”, Trans. Am. Nucl Soc., 69, 470,1995. 4. Glassner, A.S., Graphics Gems, Academic Press, Inc., San Diego, California, pp. 246-248,1990. 5. Bresenham, J.E., “Algorithm for Computer Control of a Digital Plotter”, IBM System Journal, Vol. 4, pp. 25-30,1965. 6. Foley, J.D.; van Dam, A.; Feiner, S.K.; Hughes, J.F.; Phillips, R.L., Introduction to Computer Graphics, Addison-Wesley Publishing Company, Inc., Reading, Massachusetts, pp. 70-76,1994. 7. Heckbert, P.S., Graphics Gems Volume IV, Academic Press, Inc., San Diego, California, pp. 366-368,1994. 8. Wheeler, FJ.; Wessol, D.E.; Wemple, CA. The rtt_MC Software Module for Treatment Planning. In=Venhuizen, J.R., ed. INEEL BNCT Research Program Annual Report 1996. Idaho National Engineering and Environmental Laboratory, INEL/EXT-97-00319,1997. 9. Wheeler, FJ.; Parsons, D.K.; Nigg, D.W.; Wessol, D.E.; Miller, L.G.; Fairchild, R.G., Reactor Physics Design for the Brookhaven Medical Research Reactor Epithermal Neutron Source. In: Harling, O.K., Bernard, J.A., Zamenhof, R., eds., Neutron Beam Design, Development, and Performance for Neutron Capture Therapy. New York: Plenum Press, pp. 83-95,1990. 10. E. Straker, P. Stevens, D. Irving, and V. Cain, “MORSE-CG, General Purpose Monte-Carlo Multigroup Neutron and Gamma-Ray Transport Code with Combinatorial Geometry”, Radiation Shielding Information Center fORNLl Report Number CCC-203.1976. http://id.inel.gov/cart/rtpe-manual/secOO.html http://id.inel.gov/cart/rtpe-manual/secOO.html 49 APPENDIX Three-Dimensional Algorithm to Examine Representative Univels along Line /* This algorithm w ill terminate a f te r looping through NUMSTEPS univels. * The code given is for the case of x varying most rapidly, positive * i n i t i a l coordinates, and positive d irec tion vector components, Other * cases are sim ilar but would c lu t te r the logic. * / void line3DFIR(double * p ts , double * d ir , in t NUMSTEPS) { in t i , incrx l, incrx2, incrxy, incrxz, px, py, pz, errxy, e rrxz ; double x, y, z, xd ir, ydir, zd ir, mxy, mxz, fcx, fy, fz; x = *pts++; y = *pts++; z = *pts; xdir = dir++; ydir = *dir++; zd ir = * d ir ; mxy = y d i r /x d i r ; mxz = zd ir /xd ir ; incrxl = mxy*LARGE_INT; incrxy = incrxl-LARGE_INT; incrx2 = mxz*LARGE_INT; incrxz = in crx2-L ARGE_INT; /* Position (x, y, z) */ /* Direction (xdir, ydir, zd ir V /* Slopes */ /* Increments for e rror terms V px = INT_FLOOR_FCN(x+0,5); /* Center ray for primary direction */ fcx = (double)px + 0.5; fy = y + mxy*( f cx-x); fz = z + mxz*( fcx-x); py = INT_FLOOR_FCN( fy ); pz = INT_FLOOR_FCN( f z ); /* Compute i n i t i a l e rror terms */ errxy = ( fy - (double)py)*LARGE_INT+incrxy; errxz = ( fz - (double)pz)*LARGE_INT+incrxz; /* Examine the material at the i n i t i a l position */ EXAMINE_3D (INT_FLOOR_FCN (x ), INT_FLOOR_FCN (y ), INT_.FLOOR_FCN (z )); /* Actual tracking loop begins h e r e ------------------------------------*/ for (i=0; i=0) { /* Adjust y and xy erro rs */ errxy+=incrxy; 50 py++; } else { errxy+=incrxl; } i f (errxz>=0) { /* Adjust z and xz e rro rs */ errxz+=incrxz; pz++; } else { errxz+=incrx2; } } } 51 Three-Dimensional Algorithm to Examine All Univels Intersected by Line /* This algorithm will terminate a f te r looping through NUMSTEPS univels. * The code is for the general case of an a rb i tra ry s t a r t position and * an a rb i tra ry unit vector d i re c t io n . * Governing Equations: x*dy - y*dx + c = 0 */ * x*dz - z*dx + d = 0 */ * y*dz - z*dy + e = 0 */ void line3DFIA(double * pts, double * d ir , in t NUMSTEPS) { double x, y, z, xdir, ydir, zd ir, xposdir, yposdir, zposdir, to_xy_floor, to_xz_floor, to_yz_floor, to_xy, to_xz, to_yz, inv_xdir, inv_ydir, inv_zdir; in t j , px, py, pz, incrx, incry, incrz, errxy, e rrx z , e rry z , move_x, move_y, move_z; x = pts++; /* Position (x, y, z) V y = *pts++; z = *pts; xdir = xposdir = dir++; /* Direction (xdir, ydir, zd ir */ ydir = yposdir = *dir++; zdir = zposdir = * d ir ; i f (xposdir<0,0) xposdir = -xposdir; i f (yposdir<0.0) yposdir = -yposdir; i f (zposdir<0.0) zposdir = -zposdir; incrx = xposdir*LARGE_INT; /* Increments for error terms */ incry = yposdir*LARGE_INT; incrz = zposdir*LARGE_INT; px = INT_FLOOR_FCN(x ); /* Current array position V py = INT_FLOOR_FCN(y); pz = INT_FLOOR_FCN(z ); to_yz_fIoor = x-( double )px; /•* Distances from edges V to_xz_floor = y - (double)py; to_xy_fIoor = z-(double)pz; i f (xdir>0.0) { /* I n i t i a l i z e for x V to_yz = I . 0-to_yz_floor; move_x = I; inv_xdir = I , 0/xdir; } else i f (xdir<0.0) { to_yz = to_yz_floor; move_x = - I ; inv_xdir = I . 0/xdir; } e lse { /* uncommon special case V to_yz = to_yz_floor; i f (to_yz<0.5) to_yz = I .0-to_yz; 52 move_x = 0; inv_xdir = 0.0; } i f (ydir>0,0) { /* I n i t i a l i z e for y */ to_xz = I . 0-to_xz_floor; move_y = I; inv_ydir = I .0 /y d ir ; } else i f (ydir<0.0) { to_xz = to_xz_floor; move_y = - I ; inv_ydir = I . 0/ydir; } e lse { /* uncommon special case V to_xz = to_xz_floor; i f ( to_xz<0.5) to_xz = I . 0-to_xz; move_y = 0; inv_ydir = 0.0; } i f (zdir>0.0) { /* I n i t i a l i z e for z V to_xy = I .0-to_xy_floor; move_z = I; inv jzd ir = I .0 /zdir; } else i f (zdir<0,0) { to_xy = to_xy_floor; move_z = - I ; inv_zdir = I , 0/zdir; } else { /* uncommon special case */ to_xy = to_xy_floor; i f (to_xy<0.5) to_xy = I .0-to_xy; to_xy = 1.0; move_z = 0; inv jzd ir = 0.0; } /* Compute i n i t i a l e rro r terms */ errxy = (yposdir*to_yz - xposdir*to_xz)*LARGE_INT; errxz = (zposdir*to_yz - xposdir*tp_xy)*LARGE_INT; erryz = (zposdir*to_xz - yposdir*to_xy)*LARGE_INT; /* Examine the material at the i n i t i a l position V EXAMINE_3D(px, py, pz); /* Actual tracking loop begins here ------------------------------------V for (i=0; i0) { /* don 't change y V i f (errxz>0) { /* don 't change y or x */ pz+=move_z; errxz-=incrx; erryz-=incry; } else { /* don 't change y or z V 53 px+=move_x; errxy+=incry; errxz+=incrz; } } else { i f (errxy>0) { /? don11 change z */ / * don't change z or x V py+=move_y; errxy-=incrx; erryz+=incrz; } e lse { px+=move_x; errxy+=incry; errxz+=incrz; } /* don't change z or y */ } /* Examine current m aterial; typ ica lly stop when d iffe ren t */ EXAMINE_3D(px, py, pz); } } 54 Timing Comparisons of Various Algorithms through Same Nmnher of Elements The following table compares various tracking algorithms with varying numbers of iterations based on the number of elements (pixels or univels) examined. The number of elements examined varies from I to 1000. To give representative values, each time measured is the sum of the time required to track 100,000 rays through the given number of elements. Times given are in seconds. - Elems. 2 D 2D 2D 2D 2D 2D 3D 3D 3D 3D 3D 3D DR FIR FFR RA FIA FFA DR FIR FFR RA FIA FFA I 0.05 0 .36 0 .27 0 .04 0.26 0.15 0 .07 0.55 0 .37 0 .07 . 0.43 0 .26 2 0 .07 0.37 0 .32 0 .06 0.28 0.18 0 .09 0.57 0.48 0.09 0.45 0.31 3 0 .08 0.38 0.38 0.08 0.29 0 .20 0.11 0.59 0 .59 0.11 0.48 0.37 4 0 .10 0.39 0.45 0 .10 0.31 0.23 0 .12 0.61 0.68 0 .14 0.50 0.43 5 0.11 0.41 0 .50 0.11 0.32 0 .26 0 .16 0.62 0.78 0 .17 0 .54 0.50 6 0.13 0.43 0.58 0 .12 0.35 0 .28 0 .17 0.66 0:91 0 .19 0.55 0 .57 7 0 .14 0.44 0 .64 0 .14 0.35 0.31 0 .19 0.67 1.01 0 .22 0.58 0 .62 8 0.15 0.45 0 .70 0.15 0.37 0.33 0.21 0.68 1.11 0 .24 0.59 0.68 9 0 .17 0.47 0 .76 0 .17 0.38 0.35 0.23 0.70 1.22 0 .26 0.62 0.73 10 0 .19 0.48 0 .82 0.18 0.40 0 .37 0 .24 0.72 1.32 0 .29 0.63 0.79 20 0 .30 0.59 1.42 0 .30 0.52 0 .59 0 .40 0.87 2.35 0 .50 0.84 1.36 30 0 .42 0.70 2.01 0.43 0.65 0 .80 0.55 1.01 3 .37 0.71 1.03 1.93 40 0 .54 0.80 2.61 0 .56 0.76 1.01 0 .70 1.15 4 .39 0.91 1.23 2 .49 50 0 .66 0.91 3 .20 0 .68 0.88 1.21 0 .84 1.29 5.42 1.11 1.42 3.05 60 0.78 1.01 3 .80 0 .80 1.00 1.42 0 .99 1.42 6 .44 1.32 1.60 3 .62 70 0 .89 1.11 4 .40 0 .92 1.12 1.63 1.14 1.56 7.45 1.51 1.79 4 .1 8 80 1.01 1.2 2 - 4.99 1.04 1.24 1.84 1.28 1.71 8.48 1.71 1.97 4 .74 90 1.12 1.32 5.58 1.17 1.36 2 .04 1.43 1.84 9.51 1.91 2 .16 5.30 100 1.24 1.43 6.18 1.29 1.47 2 .24 1.58 1.97 10.53 2.11 2.35 5.86 200 2 .39 2.48 12.13 2.51 2.67 4.31 3.03 3.35 20 .76 4 .09 4 .19 11.45 300 3 .54 3.52 18.08 3.73 3 .84 6 .36 4 .48 4 .72 31 .00 6 .06 6.03 17.05 400 4 .70 4 .57 24 .04 4.95 5.03 8.42 5.93 6.08 41 .24 8.03 7.86 22 .64 500 5.85 5.62 30 .00 6 .17 6.22 10.47 7 .38 7.45 51.47 9 .99 9.71 28.23 600 7 .00 6.66 35.95 7 .39 7.39 12.52 8.83 8.82 61 .70 11.96 11.53 33 .82 700 8.15 7.71 41 .90 8.61 8.58 14.58 10.28 10.19 71 .94 13.93 13.37 39.40 800 9.31 8.76 47 .85 9.83 9 .77 16.64 11.73 11.55 82.18 15.89 15.21 44.98 900 10.46 9.80 53.81 11.05 10.95 18.69 13.18 12.91 92.42 17.85 17.04 50.58 1000 11.62 10.85 59.76 12.27 12.13 20.75 14.63 14.28 102.60 19.83 18.87 56.17