{"id":175,"date":"2018-02-10T19:01:02","date_gmt":"2018-02-10T18:01:02","guid":{"rendered":"http:\/\/wp.unil.ch\/gaia\/?page_id=175"},"modified":"2019-12-10T14:24:04","modified_gmt":"2019-12-10T13:24:04","slug":"ds-matlab","status":"publish","type":"page","link":"https:\/\/wp.unil.ch\/gaia\/mps\/ds-matlab\/","title":{"rendered":"A Matlab code for DS"},"content":{"rendered":"

Overview<\/strong><\/h2>\n

The Matlab\u00ae<\/sup> code presented below is a very simplified and schematic implementation of the Direct Sampling (DS) algorithm (more information on the algorithm here<\/a>). 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\u2019s 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).<\/p>\n

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.<\/p>\n

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.<\/p>\n

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.<\/p>\n

The Matlab<\/strong>\u00ae<\/sup><\/strong> DS code<\/strong><\/h2>\n
%DO NOT USE FOR BENCHMARKING\r\n\r\n%% Initialization\r\n clear; home;\r\n %Performs a multiple-points simulation by Direct Sampling.\r\n %The training image is generated using objects.\r\n %Parameters are:\r\n simul_size = [80 80];\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %size of the simulation: y x\r\n ti_size = [100 100];\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %size of the ti: y x\r\n template = [9 9];\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %size of the template: y x\r\n fract_of_ti_to_scan = 0.1; %maximum fraction of the ti to scan\r\n thr = 0.0;\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %threshold position (between 0 and 1)\r\n \r\n %% Generation of a training image and display\r\n ti=zeros(ti_size(1),ti_size(2)); [xx,yy]=meshgrid(1:ti_size(1),1:ti_size(2));\r\n for i=1:35\r\n \u00a0\u00a0 obj_param=[random('unif',1,100) random('unif',1,100),random('unif',4,8)];\r\n \u00a0\u00a0 d=sqrt((xx-obj_param(1)).^2 + (yy-obj_param(2)).^2);\r\n \u00a0\u00a0 ti(d <= obj_param(3))=1;\r\n end\r\n \r\n figure(1); clf; subplot(1,2,1); colormap gray;\r\n imagesc(ti); title('training image'); axis equal tight xy; drawnow;\r\n H=gca; set(H,'position',[0, 0.25, 0.5 0.5]); tic;\r\n \r\n %% DS Simulation\r\n %defining shifts related to the template size\r\n yshift = floor(template(1)\/2); xshift = floor(template(2)\/2);\r\n %reducing the size of the ti to avoid scanning outside of it\r\n ti_size(1) = size(ti,1)-(2*yshift); ti_size(2) = size(ti,2)-(2*xshift);\r\n %creating empty simulation with a wrapping of NaNs of the size \"shift\"\r\n simul = nan(simul_size(1)+(2*yshift),simul_size(2)+(2*xshift));\r\n %defining path in ti and in simulation\r\n path_ti = randperm(ti_size(1)*ti_size(2));\r\n path_sim = randperm(simul_size(1)*simul_size(2));\r\n sizesim = size(path_sim,2);\r\n progress = 0; tinod = 0;\r\n \r\n %looping simulation nodes\r\n for simnod = 1:sizesim\r\n \r\n \u00a0\u00a0 progress_current=ceil((simnod*100)\/sizesim);\r\n \u00a0\u00a0 if progress_current>progress\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 progress=progress_current; disp([num2str(progress),' % completed'])\r\n \u00a0\u00a0 end\r\n \u00a0\u00a0 \r\n \u00a0\u00a0 %find node in the simulation grid\r\n \u00a0\u00a0 xsim = ceil(path_sim(simnod)\/simul_size(1));\r\n \u00a0\u00a0 ysim = path_sim(simnod)-((xsim-1)*simul_size(1));\r\n \r\n \u00a0\u00a0 %define the point and shifting it to avoid scanning NaNs\r\n \u00a0\u00a0 point_sim = [ysim+yshift xsim+xshift];\r\n \r\n \u00a0\u00a0 %define data event at simulated point\r\n \u00a0\u00a0 data_event_sim = simul(point_sim(1)-yshift:point_sim(1)+yshift ,...\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 point_sim(2)-xshift:point_sim(2)+xshift);\r\n \r\n \u00a0\u00a0 %scan the ti\r\n \u00a0\u00a0 mindist = inf; %initial best distance\r\n \u00a0\u00a0 tries = 0;\u00a0\u00a0\u00a0\u00a0 %counter of attempts\r\n \u00a0\u00a0 max_scan = size(path_ti,2)*fract_of_ti_to_scan; %max number of scans\r\n \r\n \u00a0\u00a0 %reducing the data event to its informed nodes\r\n \u00a0\u00a0 no_data_indicator = isfinite(data_event_sim);\r\n \u00a0\u00a0 data_event_sim = data_event_sim(no_data_indicator); \r\n \r\n \u00a0\u00a0 while 1==1 %scan the ti\r\n \r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 tinod = tinod+1; tries = tries+1;\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %if arriving at the end of the path, restart from the beginning\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if (tinod > size(path_ti,2))\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 tinod = 1;\r\n \u00a0\u00a0\u00a0\u00a0\u00a0end\r\n \r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %find the point in the ti\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 xti = ceil(path_ti(tinod)\/ti_size(1));\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 yti = path_ti(tinod)-((xti-1)*ti_size(1));\r\n \r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %find scanned point and data event\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 point_ti = [yti+yshift xti+xshift];\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 data_event_ti = ti(point_ti(1)-yshift:point_ti(1)+yshift ,...\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 point_ti(2)-xshift:point_ti(2)+xshift);\r\n \r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %if template is totally unknown, take the first value\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if sum(no_data_indicator(:)) == 0;\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 simul(point_sim(1),point_sim(2)) = ti(point_ti(1),point_ti(2));\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 break\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 end\r\n \r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %find the data event at this point in the ti\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 data_event_ti = data_event_ti(no_data_indicator);\r\n \r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %evaluate the distance between both data events\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %using another distance here allows using continuous variable\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 distance=mean(data_event_sim~=data_event_ti);\r\n \r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %if distance under threshold, the point is accepted\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if distance <= thr\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 simul(point_sim(1),point_sim(2)) = ti(point_ti(1),point_ti(2));\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 break\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 else\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %check if the distance is under the minimum found so far\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if distance < mindist\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 mindist = distance;\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 bestpoint = point_ti;\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 end\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 %if max_scan nodes have been scanned, take the best point found\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 if tries > max_scan\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 simul(point_sim(1),point_sim(2))=ti(bestpoint(1),bestpoint(2));\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 break\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 end\r\n \u00a0\u00a0\u00a0\u00a0\u00a0\u00a0 end\r\n \u00a0\u00a0 end\r\n end\r\n \r\n % remove the \"wrapping\" of NaNs\r\n toc; simul = simul(yshift+1:end-yshift,xshift+1:end-xshift);\r\n \r\n %show simul\r\n subplot(1,2,2); imagesc(simul); title('simulation'); axis equal tight xy; \r\n H=gca; set(H,'position',[0.5, 0.25, 0.5*simul_size(1)\/ti_size(1) ...\r\n \u00a0\u00a0 0.5*simul_size(2)\/ti_size(2)])\r\n \r\n<\/pre>\n
Using the code<\/strong><\/h2>\n
For convenience, the lines above is available in the file DS_demo.m<\/a>. 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.<\/p>\n
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).<\/p>\n
Outputs<\/strong><\/h2>\n
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.<\/p>\n
\"\"<\/p>\n
Figure <\/strong>1<\/strong> Training image and simulation produced by the Matlab DS code.<\/p>\n","protected":false},"excerpt":{"rendered":"
Overview The Matlab\u00ae 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 … <\/p>\n
Continue reading “A Matlab code for DS”<\/span><\/a><\/p>\n","protected":false},"author":1001689,"featured_media":0,"parent":197,"menu_order":3,"comment_status":"closed","ping_status":"closed","template":"page-templates\/full-width-page.php","meta":{"_seopress_robots_primary_cat":"","_seopress_titles_title":"","_seopress_titles_desc":"","_seopress_robots_index":"","footnotes":""},"class_list":["post-175","page","type-page","status-publish"],"_links":{"self":[{"href":"https:\/\/wp.unil.ch\/gaia\/wp-json\/wp\/v2\/pages\/175","targetHints":{"allow":["GET"]}}],"collection":[{"href":"https:\/\/wp.unil.ch\/gaia\/wp-json\/wp\/v2\/pages"}],"about":[{"href":"https:\/\/wp.unil.ch\/gaia\/wp-json\/wp\/v2\/types\/page"}],"author":[{"embeddable":true,"href":"https:\/\/wp.unil.ch\/gaia\/wp-json\/wp\/v2\/users\/1001689"}],"replies":[{"embeddable":true,"href":"https:\/\/wp.unil.ch\/gaia\/wp-json\/wp\/v2\/comments?post=175"}],"version-history":[{"count":0,"href":"https:\/\/wp.unil.ch\/gaia\/wp-json\/wp\/v2\/pages\/175\/revisions"}],"up":[{"embeddable":true,"href":"https:\/\/wp.unil.ch\/gaia\/wp-json\/wp\/v2\/pages\/197"}],"wp:attachment":[{"href":"https:\/\/wp.unil.ch\/gaia\/wp-json\/wp\/v2\/media?parent=175"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}