Fast Variable Density Poisson-Disc Sample Generation with Directional Variation
Nicholas Dwork, Corey A. Baron, Ethan M. I. Johnson, Daniel O'Connor, John M. Pauly, Peder E. Z. Larson
FF AST V ARIABLE D ENSITY P OISSON -D ISC S AMPLE G ENERATION WITH D IRECTIONAL V ARIATION
A P
REPRINT
Nicholas Dwork ∗ Department of Radiology and Biomedical ImagingUniversity of California in San Francisco
Corey A. Baron
Robarts Research InstituteThe University of Western Ontario
Ethan M. I. Johnson
Department of Biomedical EngineeringNorthwestern University
Daniel O’Connor
Department of Mathematics and StatisticsUniversity of San Francisco
John M. Pauly
Department of Electrical EngineeringStanford University
Peder E. Z. Larson
Department of Radiology and Biomedical ImagingUniversity of California in San FranciscoApril 16, 2020 A BSTRACT
We present a fast method for generating random samples according to a variable density Poisson-discdistribution. A minimum threshold distance is used to create a background grid array for keepingtrack of those points that might affect any new candidate point; this reduces the number of conflictsthat must be checked before acceptance of a new point, thus reducing the number of computationsrequired. We demonstrate the algorithm’s ability to generate variable density Poisson-disc samplingpatterns according to a parameterized function, including patterns where the variations in density area function of direction. We further show that these sampling patterns are appropriate for compressedsensing applications. Finally, we present a method to generate patterns with a specific accelerationrate. K eywords Poisson Disc · Parallel Imaging · Compressed Sensing
In MRI, multiple coils and compressed sensing have both reduced the number of samples required to generate di-agnostic quality imagery. The multiple coils provide additional spatial encoding that is used to interpolate missingk-space data points (a technique commonly called parallel imaging ) [1, 2, 3]. Compressed sensing takes advantage ofthe a priori knowledge that most of the values of the image are approximately after a sparsifying linear transforma-tion (e.g. a Daubechies Wavelet transform). When the system matrix satisfies specific properties (e.g. the RestrictedIsometry Principal, the Restricted Isometry Principal in Levels, or the Mutual Coherence Conditions) then the erroron the final image is bounded [4, 5, 6]. Remarkably, these conditions can often be achieved with a random samplingpattern [7]. Compressed sensing has been used in MRI with great success [8, 9, 10].When combining multi-coil imaging with compressed sensing, one wants to employ a random sampling pattern tosatisfy the compressed sensing requirements, but still keep samples far enough from each other to take advantageof the spatial encoding of the multiple coils. A poisson-disc sampling pattern can be used to simultaneously satisfythese properties [9]. It is desirable that the pattern generation algorithm be fast in order to permit investigation ofdifferent sampling distributions and determine the advantages and disadvantages of each one. Furthermore, it should ∗ a r X i v : . [ ee ss . I V ] A p r PREPRINT - A
PRIL
16, 2020accommodate densities that depend on direction to account for different coil configurations. And finally, the samplingpattern should satisfy a desired overall acceleration factor.The simplest dart throwing algorithm for generating a poisson disc sampling pattern (randomly choose a point, verifythat the point is not too close to any existing point, repeat) is a slow process [11]. Other methods are more efficient, butimpose additional requirements (e.g. a mesh defining a surface of interest, or a tiling of the space where the densityalong the edges of the tiles may be noticeably different) [12]. These are confounding effects that are not requiredfor sample generation in MRI, where the region of k-space of interest is a simple rectangular subset of a Euclideanspace. In [13], Bridson described a fast O ( n ) algorithm for generating a poisson disc sampling pattern with a constantdensity. In [14], Tulleken adapted this algorithm to accommodate a variable density sampling pattern based on a priori knowledge of the largest poisson-disc radius parameter. Tulleken’s method is O (cid:0) n (cid:1) . Notably, this method cannotaccommodate a sampling density that depends on direction.In many cases of sample generation with MRI, we know the minimum distance between samples a priori . In thispaper, we alter the methods of [13, 14] to take advantage of this knowledge. We present a faster method of generatingsamples according to a variable density Poisson-disc distribution in a rectangular subset of a Euclidean space withan arbitrary number of dimensions. We make three novel contributions for generating variable density poisson discsampling patterns to be used in MRI. • We present a more computationally efficient and faster method than the state of the art. • We present a computationally efficient method to accommodate rotationally asymmetric sampling density,permitting more acceleration in one direction than another. • We present an automatic method to generate a sampling pattern with a desired overall acceleration rate, evenwith an asymmetric sampling density.
The method of [13] reduces the computational time dramatically over the dart throwing algorithm by utilizing abackground grid. Assuming a fixed Poisson-disc parameter r , the method partitions the space into a set of cubeswhere the edges have length r/ √ d (where d is the number of dimensions of the space). The cube edge is set to thislength so that the cube’s diagonal has length r ; thus, each cell can contain at most point. For each new point, then,one need not check the distance to all other points. Instead, only those points that are indexed in grid cells within range(as illustrated in the Fig. 1a, where r = r max ) need to be checked [14].Figure 1: a) For a new point (shown with a black dot), only those cells intersected by the red circles (shaded blue)can contain points close enough to violate the distance threshold. If the point were in n dimensions instead of dimensions, the length of each side of all cubes would be r max / √ n . b) A new candidate point y j with Poisson-discparameter r j is illustrated with the black dot. In this figure, ceiling ( r j /r min ) = 4 . The red circles have diameter r min and are tangent to the corner of the cell containing y j . The blue squares are those cells that could contain pointsthat are within r j of y j . If the point were in n dimensions instead of dimensions, the length of each side of all cubeswould be r min / √ n .The set of valid points is initialized by choosing a point x at random in the space. A list of point indices (called theactive list) is initialized with index . An array of size equal to the background grid is created and initialized so that all2 PREPRINT - A
PRIL
16, 2020values are (meaning that the corresponding cell does not contain any points). The element of the background arraythat corresponds to the cell containing x is set to . Computation then proceeds as detailed in Alg. 1. Algorithm 1:
Fast Poisson-Disc Sampling with constant r Inputs: r , active list, background array, parameter k While: active list is not emptyChoose a point randomly from the active list, x i Create k new points uniformly at random in the spherical annulusbetween radii r and r centered on x i For each: created point y j Find the indices of those points that may be closer than r to y j using the background grid as illustrated in Fig. 1a. Find thedistance between y j and the points with those indices. If the distance between y j and all existing points is greaterthan r Add y j to the list of points, add j to the active list, and addindex j to the appropriate elements of the background grid. End IfEnd For
If none of the k points were valid, remove i from the active list. End WhileOutputs:
List of pointsIn Alg. 1, one must choose points randomly on a spherical annulus. The number of points chosen is denoted by k .In [13], Bridson suggests a value of k = 30 ; we have found that the processing is faster with comparable resultsin two dimensions when k = 10 . To choose points in two dimensions, construct a vector with angle (in radians)chosen uniformly at random on [ − π, π ) and magnitude chosen uniformly at random on [ r, r ] . To sample a pointat random in n dimensions, generate a vector with n elements where the value of each element is a realization ofa normally distributed random variable. Since the normal distribution is rotationally symmetric, every direction hasequal probability density. Then, scale this vector to a magnitude chosen uniformly at random on [ r, r ] .The method by Tulleken of [14] accommodates a variable density Poisson-disc sampling pattern (meaning that theparameter r changes as a function of location) by altering the background array so that 1) each element of the gridaccepts a list of points and 2) the cell size must be computed from the maximum possible radius. The index of eachnew point is added to the list of the background array element that corresponds to the grid cell that contains the newpoint.Though this would result in a realization of the desired sampling pattern, it is inefficient, as the following thoughtexperiment illustrates. Suppose that r is small near the center of the image and increases as the distance from thecenter increases. In this case, the largest values of r would be attained at the corners of the region. Indeed, thisis the most common use case for compressed sensing MRI applications. These values could be so large that therectangular region of interest would be divided into a small number of large grid cells, meaning that many pointswould be listed within each grid cell and all of those points would need to be checked with each new additional point.The computational cost degenerates to that of the extremely slow dart throwing algorithm. In section 2.2, we explainhow to overcome this inefficiency. For the fast algorithm, it is assumed that a positive minimum bound on r exists: r min > . This is almost certainlythe case with MRI. The k-space samples need not be closer than the inverse of the field-of-view. Moreover, withcompressed sensing, the center region of k-space is often fully sampled (meaning that samples are separated by adistance equal to the inverse of the field-of-view) [15, 9]. With these applications, the density of samples do not varyunless they are located some positive distance from the origin. Thus, for any variable density scheme that reducesthe sampling density as distance from the origin increases, the size of the fully-sampled center region can be used todetermine r min .The region of interest is partitioned into a grid of cubes where the edges all have length r min / √ n . Grid elements donot contain the indices of points that fall within their boundary. Instead, each grid element contains a list of indicesof those points that might have their threshold distance violated by a point in the grid cell. When a new point is3 PREPRINT - A
PRIL
16, 2020considered as a candidate for the sample distribution, its distance is checked against all of those points with indiceslocated in its grid cell. If the candidate point is not too close to any existing points, then its index is added to all ofthose grid cells where its distance threshold reaches (as illustrated in Fig. 1b).The set of valid points is initialized by choosing a point x at random in the space. Add index to the active listand add index to all those background grid cells that might contain points within the threshold distance of point x according Fig. 1b. Then proceed according to Alg. 2. Algorithm 2:
Fast Variable Density Poisson-Disc Sampling
Inputs: r min , active list, background array, parameter k While: active list is not emptyChoose a point randomly from the sample list, x i Identify the poisson-disc parameter for this point, r ( x i ) Create k new points uniformly at random in the spherical annulus between radii r ( x i ) and r ( x i ) centered on x i For each: created point y j Compute the distance between y j and those points with indicesin the list contained in the background grid element thatcorresponds to the cell containing y j . If the distance between y j and all existing points is greaterthan r ( x j ) Add y j to the list of points, add j to the active list, and addindex j to the appropriate element of the background array. End IfEnd For
If none of the k points were valid, remove i from the active list. End WhileOutputs:
List of pointsFor the results presented in this work, we used the following function for the poisson-disc parameter r : r γ ( x ) = (cid:107) x (cid:107) + 0 . γ , (1)where (cid:107) · (cid:107) represents the L norm. For this function, r min = 0 . /γ and r max = ( (cid:107) x corner (cid:107) + 0 . /γ where x corner is any corner of the sampling domain. The rectangular region of k-space of interest is usually the [ − . , . n cube . (An exception would be a Homodynesampling pattern, where the rectangular region is a little more than half of this cube [16].) As discussed, it may bedesirable to alter the density distribution of points as a function of direction. That is, in addition to the poisson discparameter r being a function of location, it would also be a function of radial direction to take advantage of coilplacement geometry. For a general function of radial direction, this is a computationally challenging task since itrequires determining the direction between points and evaluating this function. However, in MRI, it is usually thecase that we are interested in acceleration rates that differ along the axes of the rectangular region. For example, witha birdcage coil, the spatial encoding of each coil differs significantly in the transverse plane but differs little in thelongitudinal direction. Therefore, a higher acceleration is possible in the transverse plane than in the longitudinaldirection. And example of this can be seen in Fig. 4a.Since we are interested in acceleration rates that differ along the axes of the rectangular region, a simple trick avoidsdetermining the angle between points and evaluating the parameter function. For simplicity, we will consider a twodimensional case. Suppose one wants to undersample the second dimension by an additional factor of ν (meaningthe density of samples will be greater in the first dimension than in the second). Then, one can generate a variabledensity poisson-disc sampling pattern on the region [ − . , . × [ − . /ν, . /ν ] quickly with a radially symmetricPoisson-disc parameter using Alg. 2. After the points are generated, scale the resulting pattern by ν in the seconddimension to generate a pattern of points on [ − . , . × [ − . , . . Here, set multiplication is the Cartesian cross product. PREPRINT - A
PRIL
16, 2020
The acceleration rate of a sampling pattern is the number of samples acquired divided by the number of samplesrequired for full sampling. It is often convenient to be able to specify an overall acceleration rate and attain a corre-sponding sampling pattern. In this section, we describe a method to attain this goal.In order to do so, we require a poisson-disc parameter function r γ : R d → R parameterized by γ ≥ such that r decreases monotonically with increasing γ . (Equivalently, the overall acceleration rate increases with increasing γ .)Note that (1) satisfies this property. Then, with the computational efficiency of Alg. 2, we can now specify an overallacceleration rate α and determine the desired sampling pattern with a binary search algorithm in a reasonable amountof time. (The total computational time is the time of any single iteration multiplied by the number of iterations in Alg.3.)In order to use the binary search, one must supply bounds γ min and γ max . Since γ is positive, γ min = 0 is a lowerbound. The γ that corresponds to r min would be γ max . The complete search is specified in Algorithm 3. Algorithm 3:
Binary Search to Find Pattern for Specified Acceleration Rate
Inputs:
Desired acceleration rate α , bounds γ min and γ max , tolerance tol > . Do: (cid:15) = ( γ max − γ min ) / γ mid = (cid:15) + γ min Determine the sampling pattern P using Algorithm 2 withparameter function r γ mid = α mid If:
Acceleration rate > αγ max := γ mid Else: γ min := γ mid While (cid:15) > tol
Outputs:
Pattern P Figure 2 shows variable density sampling patterns generated with Algorithm 2 using the poisson-disc parameter func-tion of (1) for γ ∈ { , , , , } and additional directional undersampling of 3. As expected, as γ increases, the overallacceleration rate (equal to the number of samples divided by the size of the domain) decreases.5 PREPRINT - A
PRIL
16, 2020Figure 2: Variable density poisson disc sampling patterns generated with the fast algorithm of Alg. 2.Table 1 compares the time required to compute the variable density poisson disc sampling pattern with Alg. 2 to thetime required by the Tulleken algorithm for the sampling patterns shown in Fig. 2. The improvement in runtime isbetween − . The algorithms were implemented in C on a 2012 Macbook Pro with a 2.5 GHz Intel i7 processor.Processing Time (ms): Alg. 2 / Tulleken γ = 50 γ = 75 γ = 100 γ = 125 γ = 150 (3,1) 10 / 20 20 / 40 40 / 70 60 / 110 110 / 150(1,1) 20 / 30 50 / 80 100 / 140 150 / 220 230 / 340(1,3) 10 / 20 20 / 40 40 / 70 70 / 110 110 / 150Table 1: Runtimes for generating variable density Poisson-disc sampling patterns with Alg. 2 and the Tullekenalgorithm. Time is reported in milliseconds. In all cases, Alg. 2 is faster by − .In Fig. 3, we present a comparison of the results from Alg. 2 to the variable density algorithm of Tulleken. For thesampling patterns with γ = 150 and without any additional directional undersampling, we created a Voronoi partitionof the domain and plotted the area of each cell versus distance of the point from the origin. Fig. 3 demonstrates thatthe distributions of area versus distance are very similar.Figures 4 and 5 show how the sampling mask could be used with MRI data with knee and ankle data, respectively. Thedata of Fig. 4 was taken from mridata.org [17]. The data for Fig. 5 was collected with a clinical 3 Tesla scannerand an 8-channel ankle coil. Both of these datasets consist of fully sampled data with two dimensions of phase-encodes and a single dimension of readout. An inverse Fast Fourier Transform was applied in the readout directionplacing the data in a hybrid space; further processing was only done on a single slice. The SAKE+L1-ESPIRiTalgorithm was used to reconstruct data that was retrospectively subsampled with the variable density poisson discsampling masks [18, 19]. Subfigure (a) shows the fully sampled reconstruction from each individual coil, (b) showsthe sum-of-squares reconstruction with fully sampled data, and (c) shows the SAKE+L1-ESPIRiT reconstruction fromretrospectively subsampled data. The SAKE+L1-ESPIRiT reconstruction is very similar to the fully sampled sum-of-squares reconstruction.Figure 6 shows patterns generated with specific accelerations using Alg. 3. Both the Tulleken and Alg. 2 wereused as the underlying method to determine the sampling pattern. A tolerance of . in the acceleration factor was6 PREPRINT - A
PRIL
16, 2020Figure 3: Area of the Voronoi cell plotted against distance from the origin for each sample point with γ = 150 forthe fast algorithm (blue) and the Tulleken algorithm (red). The distribution of points of the two algorithms is similar,indicating that they are generating sampling patterns of similar quality.Figure 4: Fully sampled reconstructions of an 8-coil extremity coil. The sensitivities of the coils in the ante-rior/posterior differ, whereas the sensitivity in the superior/inferior direction are approximately uniform for all coils.a) The fully sampled reconstructions for each of the 8 coils. b) The sum-of-squares reconstruction with fully sampleddata. c) The SAKE+L1-ESPIRiT reconstruction of the data retrospectively subsampled with the sampling pattern offigure 2 with γ = 125 and undersampling of (3,1).permitted. Indeed, sampling patterns with the specified acceleration rate were generated. And, the sampling patternsof the Tulleken algorithm and Alg. 2 are qualitatively similar.In the spirit of reproducible research, we provide a software in C with implementations of the algorithms described.The software package can be downloaded from: https://github.com/ndwork/fastVDPD_C . In this work, we presented a fast algorithm for generating a variable density Poisson-disc sampling pattern, we showhow the method can be adapted to permit a further directional undersampling, and we presented a method for gener-ating a sampling pattern with a specific acceleration factor [20]. Having a fast algorithm for creating variable density7
PREPRINT - A
PRIL
16, 2020Figure 5: Fully sampled reconstructions of an 8-coil ankle coil. a) The fully sampled reconstructions for each of the 8coils. b) The sum-of-squares reconstruction with fully sampled data. c) The SAKE+L1-ESPIRiT reconstruction of thedata retrospectively subsampled with the sampling pattern of figure 2 with γ = 100 and no directional undersampling.Figure 6: Acceleration PatternsPoisson-disc sampling patterns will empower future researchers to investigate different sampling distributions and de-termine the advantages and disadvantages of each one. Additionally, different sampling patterns can be generated fordifferent collections, which may be beneficial for some applications.When selecting candidate points around an existing point, it may be the case that the candidates are unluckily chosenso that they conflict with existing points, but other candidate points could have been created that would not have. Thus,the realized set of samples is not guaranteed to be maximal (where a maximal realization is one where the points are asdense as possible) [21]. However, the probability that the realization is maximal is typically very high (for sufficientlylarge values of k ) [22].In this paper, we have shown that a simple scaling permits a computationally efficient method for generating patternswith directional undersampling. A more general transformation can be used to create variations with more interestingpatterns. For example, an Affine transformation (or even any homeomorphism) could be used to generate samplingpatterns with variable and directional sampling patterns. The applications to MRI of such a technique are not immedi-ately obvious. 8 PREPRINT - A
PRIL
16, 2020It may be possible to further reduce the speed of sample generation by using parallelization to take advantage ofmultiple processing cores [23]. Additionally, it may be useful to adapt this algorithm to generate samples according toa variable density Poisson-disc distribution on a surface [24, 22]. We leave these prospects as possibilities for futurework.
Acknowledgements
ND and JP have been supported by the National Institute of Health’s grant number P41 EB015891 and grant numberR01 HL136965. ND would like to thank the Quantitative Biosciences Institute at UCSF and the American HeartAssociation as funding sources for this work.
References [1] Daniel K Sodickson and Warren J Manning. Simultaneous acquisition of spatial harmonics (SMASH): fastimaging with radiofrequency coil arrays.
Magnetic resonance in medicine , 38(4):591–603, 1997.[2] Klaas P Pruessmann, Markus Weiger, Markus B Scheidegger, and Peter Boesiger. SENSE: sensitivity encodingfor fast MRI.
Magnetic resonance in medicine , 42(5):952–962, 1999.[3] Mark A Griswold, Peter M Jakob, Robin M Heidemann, Mathias Nittka, Vladimir Jellus, Jianmin Wang,Berthold Kiefer, and Axel Haase. Generalized autocalibrating partially parallel acquisitions (GRAPPA).
Mag-netic Resonance in Medicine , 47(6):1202–1210, 2002.[4] David L Donoho and Michael Elad. Optimally sparse representation in general (nonorthogonal) dictionaries via (cid:96) minimization. Proceedings of the National Academy of Sciences , 100(5):2197–2202, 2003.[5] Emmanuel J Candes, Justin K Romberg, and Terence Tao. Stable signal recovery from incomplete and inaccuratemeasurements.
Communications on Pure and Applied Mathematics: A Journal Issued by the Courant Instituteof Mathematical Sciences , 59(8):1207–1223, 2006.[6] Ben Adcock, Anders C Hansen, Clarice Poon, and Bogdan Roman. Breaking the coherence barrier: A newtheory for compressed sensing. In
Forum of Mathematics, Sigma , volume 5, pages 1–84. Cambridge UniversityPress, 2017.[7] Emmanuel J Cand`es and Michael B Wakin. An introduction to compressive sampling [a sensing/samplingparadigm that goes against the common knowledge in data acquisition].
IEEE signal processing magazine ,25(2):21–30, 2008.[8] Michael Lustig, David Donoho, and John M Pauly. Sparse MRI: The application of compressed sensing for rapidMR imaging.
Magnetic Resonance in Medicine , 58(6):1182–1195, 2007.[9] SS Vasanawala, MJ Murphy, Marcus T Alley, P Lai, Kurt Keutzer, John M Pauly, and Michael Lustig. Practicalparallel imaging compressed sensing MRI: Summary of two years of experience in accelerating body MRI ofpediatric patients. In , pages1039–1043. IEEE, 2011.[10] Urs Gamper, Peter Boesiger, and Sebastian Kozerke. Compressed sensing in dynamic MRI.
Magnetic Resonancein Medicine , 59(2):365–373, 2008.[11] David Cline, Stefan Jeschke, K White, Anshuman Razdan, and Peter Wonka. Dart throwing on surfaces. In
Computer Graphics Forum , volume 28 (4), pages 1217–1226. Wiley Online Library, 2009.[12] Ares Lagae and Philip Dutr´e. A comparison of methods for generating poisson disk distributions. In
ComputerGraphics Forum , volume 27 (1), pages 114–129. Wiley Online Library, 2008.[13] Robert Bridson. Fast poisson disk sampling in arbitrary dimensions. In
SIGGRAPH sketches , page 22, 2007.[14] Herman Tulleken. Poisson disk sampling.
Dev. Mag , 21:21–25, 2008.[15] Evan Levine, Bruce Daniel, Shreyas Vasanawala, Brian Hargreaves, and Manojkumar Saranathan. 3D cartesianMRI with compressed sensing and variable view sharing using complementary poisson-disc sampling.
Magneticresonance in medicine , 77(5):1774–1785, 2017.[16] Douglas C Noll, Dwight G Nishimura, and Albert Macovski. Homodyne detection in magnetic resonance imag-ing.
IEEE transactions on medical imaging , 10(2):154–163, 1991.9
PREPRINT - A
PRIL
16, 2020[17] F Ong, S Amin, S Vasanawala, and M Lustig. Mridata.org: An open archive for sharing MRI raw data. In
Proc.Intl. Soc. Mag. Reson. Med , volume 26, page 1, 2018. .[18] Peter J Shin, Peder EZ Larson, Michael A Ohliger, Michael Elad, John M Pauly, Daniel B Vigneron, and MichaelLustig. Calibrationless parallel imaging reconstruction based on structured low-rank matrix completion.
Mag-netic resonance in medicine , 72(4):959–970, 2014.[19] Martin Uecker, Peng Lai, Mark J Murphy, Patrick Virtue, Michael Elad, John M Pauly, Shreyas S Vasanawala,and Michael Lustig. ESPIRiTan eigenvalue approach to autocalibrating parallel MRI: where SENSE meetsGRAPPA.
Magnetic resonance in medicine , 71(3):990–1001, 2014.[20] Justin P Haldar and Daeun Kim. OEDIPUS: An experiment design framework for sparsity-constrained MRI.
IEEE transactions on medical imaging , 38(7):1545–1558, 2019.[21] Cheuk Yiu Ip, M Adil Yalc¸in, David Luebke, and Amitabh Varshney. Pixelpie: maximal poisson-disk samplingwith rasterization. In
Proceedings of the 5th High-Performance Graphics Conference , pages 17–26. ACM, 2013.[22] Ethan M. I. Johnson.
Techniques and Analyses for Ultrashort Echo-time Magnetic Resonance Imaging . PhDthesis, Stanford University, 2016.[23] Li-Yi Wei. Parallel poisson disk sampling.
ACM Transactions on Graphics (TOG) , 27(3):20, 2008.[24] Xiang Ying, Shi-Qing Xin, Qian Sun, and Ying He. An intrinsic algorithm for parallel poisson disk sampling onarbitrary surfaces.