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).
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.