A Matlab code for DS

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.