Overview
The Matlab® code presented below is a very simplified and schematic implementation of the Direct Sampling (DS) algorithm (more information on the algorithm here). It corresponds to the very first version of the DS, and therefore still has characteristics inherited from traditional MP simulations. For instance, it uses fixed templates that have a box shape, which is made very easy by Matlab’s matrix-oriented language. Only the core idea of the DS is implemented, and most features described in the thesis are not present (such as continuous variable simulation, possibility to accommodate flexible data events, multivariate simulation, syn-processing, parallel computing, and so on). Multiple-grids are also not implemented (as they should because the template is fixed), and therefore large structures cannot be well reproduced. The code is presented for demonstration purposes only. The algorithm is oversimplified and it should not be used as such for real applications. Moreover, the algorithm is patented and its commercial usage is illegal without authorization (this does not concern academic research and applications).
The code is extremely short, with just over 100 lines, including comments, generation of the training image and visualization of the results. Removing comments and non-essential parts leaves a code running with only 55 lines. This illustrates the simplicity of the approach.
The first section of the Matlab code is devoted to setting the simulation parameters. The sizes of the training image, simulation grid and template, the distance threshold and the maximum fraction of the TI to scan are defined. Then comes the generation of the TI. This is done by adding in an empty grid circular objects whose location and radius are random. The DS simulation is then launched, and mainly consists of two imbricated loops: the first one on all nodes of the simulated grid (this is the sequential simulation), and the second one on all nodes of the training image (the scan of the TI). The simulation is computationally feasible because the second loop is most of the time quickly interrupted.
Another implementation, written in C++ and having all the features, was used for the examples and applications displayed in the journal papers and my my thesis. Concerning performance, the calculation times are several orders of magnitudes shorter than the Matlab version.
The Matlab® DS code
%DO NOT USE FOR BENCHMARKING
%% Initialization
clear; home;
%Performs a multiple-points simulation by Direct Sampling.
%The training image is generated using objects.
%Parameters are:
simul_size = [80 80]; %size of the simulation: y x
ti_size = [100 100]; %size of the ti: y x
template = [9 9]; %size of the template: y x
fract_of_ti_to_scan = 0.1; %maximum fraction of the ti to scan
thr = 0.0; %threshold position (between 0 and 1)
%% Generation of a training image and display
ti=zeros(ti_size(1),ti_size(2)); [xx,yy]=meshgrid(1:ti_size(1),1:ti_size(2));
for i=1:35
obj_param=[random('unif',1,100) random('unif',1,100),random('unif',4,8)];
d=sqrt((xx-obj_param(1)).^2 + (yy-obj_param(2)).^2);
ti(d <= obj_param(3))=1;
end
figure(1); clf; subplot(1,2,1); colormap gray;
imagesc(ti); title('training image'); axis equal tight xy; drawnow;
H=gca; set(H,'position',[0, 0.25, 0.5 0.5]); tic;
%% DS Simulation
%defining shifts related to the template size
yshift = floor(template(1)/2); xshift = floor(template(2)/2);
%reducing the size of the ti to avoid scanning outside of it
ti_size(1) = size(ti,1)-(2*yshift); ti_size(2) = size(ti,2)-(2*xshift);
%creating empty simulation with a wrapping of NaNs of the size "shift"
simul = nan(simul_size(1)+(2*yshift),simul_size(2)+(2*xshift));
%defining path in ti and in simulation
path_ti = randperm(ti_size(1)*ti_size(2));
path_sim = randperm(simul_size(1)*simul_size(2));
sizesim = size(path_sim,2);
progress = 0; tinod = 0;
%looping simulation nodes
for simnod = 1:sizesim
progress_current=ceil((simnod*100)/sizesim);
if progress_current>progress
progress=progress_current; disp([num2str(progress),' % completed'])
end
%find node in the simulation grid
xsim = ceil(path_sim(simnod)/simul_size(1));
ysim = path_sim(simnod)-((xsim-1)*simul_size(1));
%define the point and shifting it to avoid scanning NaNs
point_sim = [ysim+yshift xsim+xshift];
%define data event at simulated point
data_event_sim = simul(point_sim(1)-yshift:point_sim(1)+yshift ,...
point_sim(2)-xshift:point_sim(2)+xshift);
%scan the ti
mindist = inf; %initial best distance
tries = 0; %counter of attempts
max_scan = size(path_ti,2)*fract_of_ti_to_scan; %max number of scans
%reducing the data event to its informed nodes
no_data_indicator = isfinite(data_event_sim);
data_event_sim = data_event_sim(no_data_indicator);
while 1==1 %scan the ti
tinod = tinod+1; tries = tries+1;
%if arriving at the end of the path, restart from the beginning
if (tinod > size(path_ti,2))
tinod = 1;
end
%find the point in the ti
xti = ceil(path_ti(tinod)/ti_size(1));
yti = path_ti(tinod)-((xti-1)*ti_size(1));
%find scanned point and data event
point_ti = [yti+yshift xti+xshift];
data_event_ti = ti(point_ti(1)-yshift:point_ti(1)+yshift ,...
point_ti(2)-xshift:point_ti(2)+xshift);
%if template is totally unknown, take the first value
if sum(no_data_indicator(:)) == 0;
simul(point_sim(1),point_sim(2)) = ti(point_ti(1),point_ti(2));
break
end
%find the data event at this point in the ti
data_event_ti = data_event_ti(no_data_indicator);
%evaluate the distance between both data events
%using another distance here allows using continuous variable
distance=mean(data_event_sim~=data_event_ti);
%if distance under threshold, the point is accepted
if distance <= thr
simul(point_sim(1),point_sim(2)) = ti(point_ti(1),point_ti(2));
break
else
%check if the distance is under the minimum found so far
if distance < mindist
mindist = distance;
bestpoint = point_ti;
end
%if max_scan nodes have been scanned, take the best point found
if tries > max_scan
simul(point_sim(1),point_sim(2))=ti(bestpoint(1),bestpoint(2));
break
end
end
end
end
% remove the "wrapping" of NaNs
toc; simul = simul(yshift+1:end-yshift,xshift+1:end-xshift);
%show simul
subplot(1,2,2); imagesc(simul); title('simulation'); axis equal tight xy;
H=gca; set(H,'position',[0.5, 0.25, 0.5*simul_size(1)/ti_size(1) ...
0.5*simul_size(2)/ti_size(2)])
Using the code
For convenience, the lines above is available in the file DS_demo.m. The code can simply be pasted in the Matlab command window and it will execute. Alternatively, it can be copied in the editor window, where modifications on the parameters can take place. For example, increasing the scanned fraction of the TI dramatically increases simulation time. Increasing the threshold parameter reduces CPU cost, but also decreases simulations quality.
Changing the density and shape of the objects can be done with slight alterations of the code. Using a continuous (for example multiGaussian) TI would necessitate changing the measure of distance when scanning the TI, and adopting another distances (as proposed in theory).
Outputs
Running the code as presented above takes about a minute and produces an image similar to Figure 1, with the generated TI on the left and one unconditional multiple-points simulation on the right.
Figure 1 Training image and simulation produced by the Matlab DS code.