bpdq_generate_angiogram : Generate artificial Angiogram signal
function result = bpdq_generate_angiogram(m,n,num,seed)
Inputs:
m,n - size of desired synthetic angiogram
num - number of ellipses
seed - random seed
Outputs:
result - computed synthetic angiogram
Adapted from code from Ewout van den Berg and Michael P. Friedlander
GENERATEANGIOGRAM Generate artificial Angiogram signal.
GENERATEANGIOGRAM(M,N) generates an M by N image filled with 18
non-overlapping ellipses.
Based on the SparseMRI implementation by M. Lustig, 2007.
References
M. Lustig, D.L. Donoho, and J.M. Pauly, Sparse MRI: The
application of compressed sensing for rapid MR imaging,
Submitted to Magnetic Resonance in Medicine, 2007.
This file is part of Sparco Toolbox (http://www.cs.ubc.ca/labs/scl/sparco/)
Copyright (C) 2008, Ewout van den Berg and Michael P. Friedlander
(See the notice at the end of the file.)0001 function result = bpdq_generate_angiogram(m,n,num,seed) 0002 % bpdq_generate_angiogram : Generate artificial Angiogram signal 0003 % 0004 % function result = bpdq_generate_angiogram(m,n,num,seed) 0005 % 0006 % Inputs: 0007 % m,n - size of desired synthetic angiogram 0008 % num - number of ellipses 0009 % seed - random seed 0010 % Outputs: 0011 % result - computed synthetic angiogram 0012 % 0013 % Adapted from code from Ewout van den Berg and Michael P. Friedlander 0014 % 0015 % GENERATEANGIOGRAM Generate artificial Angiogram signal. 0016 % 0017 % GENERATEANGIOGRAM(M,N) generates an M by N image filled with 18 0018 % non-overlapping ellipses. 0019 % 0020 % Based on the SparseMRI implementation by M. Lustig, 2007. 0021 % 0022 % References 0023 % M. Lustig, D.L. Donoho, and J.M. Pauly, Sparse MRI: The 0024 % application of compressed sensing for rapid MR imaging, 0025 % Submitted to Magnetic Resonance in Medicine, 2007. 0026 % 0027 % This file is part of Sparco Toolbox (http://www.cs.ubc.ca/labs/scl/sparco/) 0028 % Copyright (C) 2008, Ewout van den Berg and Michael P. Friedlander 0029 % (See the notice at the end of the file.) 0030 0031 0032 % Initialise the random number generator 0033 if~exist('seed','var') 0034 seed=16000; 0035 end 0036 0037 rand('twister',seed); 0038 0039 % Settings 0040 imgSize = [m,n]; 0041 nMagnitude = 3; 0042 %sizes = [1,3,6]; 0043 sizes = [1,3,6]; 0044 magnitudes = (1:nMagnitude) / nMagnitude; 0045 0046 % Initialise the vectors 0047 n = round(length(sizes) * (length(sizes) + 1) / 2); 0048 sx = []; sy = []; 0049 mg = zeros(n * nMagnitude,1); 0050 0051 % Set up the x- and y-radii and intensity values 0052 mg = repmat(magnitudes,n,1); 0053 for i=1:length(sizes) 0054 sx = [sx; repmat(sizes(i), i,nMagnitude)]; 0055 sy = [sy; repmat(sizes(1:i)',1,nMagnitude)]; 0056 end 0057 sx = sx(:) * (imgSize(2) / 100.0); 0058 sy = sy(:) * (imgSize(1) / 100.0); 0059 mg = mg(:); 0060 0061 sx=sx(1:num); 0062 sy=sy(1:num); 0063 0064 % --------------------------------------------------------------------- 0065 % Generate set of non-overlapping ellipsoids. 0066 % Allow at most 'maxfail' occurences of overlapping ellipsoids 0067 % to avoid potential infinite looping in dense settings. 0068 % --------------------------------------------------------------------- 0069 0070 maxfail = 1000; 0071 m = imgSize(1); 0072 n = imgSize(2); 0073 img = zeros(m,n); % Image 0074 ell = zeros(m,n); % Single ellipsoid 0075 [x,y] = meshgrid((0:n-1)-round(n/2), (0:m-1)-round(m/2)); 0076 0077 for i=1:length(sx) 0078 angle = 2*pi*(rand(1)-0.5); 0079 while (maxfail > 0) 0080 ell = ellipsoid_intrnl(sx(i),sy(i),angle); 0081 if sum(any(ell .* img)) == 0, break; end; % No overlap 0082 maxfail = maxfail - 1; 0083 end 0084 img = img + mg(i) * ell; 0085 end 0086 0087 result = img; 0088 0089 0090 function res = ellipsoid_intrnl(xRadius,yRadius,angle) 0091 0092 b = sqrt(2 / (1 + xRadius^2 / yRadius^2)); 0093 a = b * xRadius / yRadius; 0094 xr = a * (x * cos(angle) - y * sin(angle)); 0095 yr = b * (y * cos(angle) + x * sin(angle)); 0096 mx = max(abs(xr(:))); 0097 my = max(abs(yr(:))); 0098 0099 res = 0; 0100 while (sum(res(:)) == 0) % While empty plot 0101 xc = 1.1 * (rand(1,1)-0.5) * mx; 0102 yc = 1.1 * (rand(1,1)-0.5) * my; 0103 res = (((xr-xc).^2 + (yr-yc).^2) < yRadius^2); 0104 end 0105 end 0106 0107 end % program 0108 0109 % The Sparco Toolbox is free software: you can redistribute it and/or modify 0110 % it under the terms of the GNU General Public License as published by 0111 % the Free Software Foundation, either version 3 of the License, or 0112 % (at your option) any later version. 0113 % 0114 % The Sparco Toolbox is distributed in the hope that it will be useful, 0115 % but WITHOUT ANY WARRANTY; without even the implied warranty of 0116 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 0117 % GNU General Public License for more details. 0118 % 0119 % You should have received a copy of the GNU General Public License 0120 % along with Sparco Toolbox. If not, see <http://www.gnu.org/licenses/>.