Displaying the output

Output files

Recall that several output (text) files may be available:

* Curves: coordinates of front points. This is a text file whose lines are in format: X Y. There is no delimiter between two fronts, even if they are associated with two different time steps.
* CurveLengths: contains, for each time step (at the end of which the front is saved) the number of points on the front (as saved in Curves).
X and Y: coordinates of mesh points (orthogonal mesh).
* Time: time discretization.
** Phi: the level set (on the whole mesh) at the last iteration (that is saved), useful to get arrival times in the fast marching method.
F: speed rates on the whole mesh at the last iteration (that is saved).
Points, Edges and Triangles: empty files that would store mesh data in case of the use of unstructured meshes (not implemented for the moment).

* Relevant for the narrow band level set method.

** Relevant for the fast marching method.

MultivacDisplay

MultivacDisplay is a small package to display results from Multivac. It is available on Multivac homepage: http://spacetown.free.fr/fronts/. Download it and expand it (tar zxvf MultivacDisplay.tgz). Notice that you need Gnuplot 3.7 or 3.8 (which requires MultivacDisplay >= 2.1) and Python 2.3 (it should work with version 2.2). All scripts are written in Python and they call Gnuplot to generate images. To generate animations (with multivac_anim), you must have netpbm (ppmtogif) and gifsicle installed too (they are in Debian packages).

Your have to compile multivac_format.cpp which is a C++ code called by the scripts. A makefile is provided: launch make and it should compile multivac_format successfully if you have a standard installation (i.e. with GNU g++). Remember that multivac_format must be available to the scripts. It means that it must be in your $PATH directories (which usually include the current directory).

In your output directory (where output files are saved), you can launch the scripts to get the images. Let's try with an example provided with Multivac: track.cpp. Launch track (as previously explained in this user's guide). Move to your output directory. You are going to generate a png file.

Type multivac_png to get help. You need to get the value of "\sqrt{delta x^2 + delta y^2}". In track.cpp, you may find that "delta x = delta y = 3.0 / 300. = 0.01". So, "\sqrt{delta x^2 + delta y^2} = 0.01414". Therefore, you should launch:

multivac_png 0.01414 1 out.png

The png file "out.png" should be displayed by any browser or image viewer. You should get the following image:

To generate an animation:

multivac_anim 0.01414 1 out.gif

You get:

Displaying with Matlab

Warning: this section is not up-to-date and Matlab scripts are not maintained anymore. It is advocated that you use MultivacDisplay instead.

All the curves saved (at different times) are written in the same file "Curves" in one stream, that is to say in one block: there is no separator between two different curves and it is impossible to tell where a curve starts and where it ends. That's why there is another file called "CurveLengths" which is created along with "Curves" and which stores the number of points for each curve, so that we can sort the points of "Curves" curve by curve and then display them. An easy way to do it is to use Matlab.

Example: displays all the saved curves.

We can then write such matlab functions to display the results of Multivac which are all saved in files described in the previous sections.

There are just two things to keep in mind. First, the different curves saved are written in the same file so we need to use the Curvelength file as seen previously. The second important thing is that in the file "Curves", a curve is a succession of points. So if a curve is in several different parts (for example 2 circles), nothing indicates how to build the curve right. So if we just ask Matlab to create segments between each set of two points, then the display will join points that belong to different parts of the curve. The solution is to test for each set of two points the size of the segment. If this size is greater than the diagonal of one cell of our mesh (width Delta_x and height Delta_y), then it means that these two points can't belong to the same part a the curve. Indeed, we know that every points belong to the mesh and that the diagonal is the greatest distance between two points of a rectangle, so two points belonging to a same part will belong to the same rectangle and will then have a distance lower than the diagonal which is sqrt( Delta_x^2 + Delta_y^2 ). These really important features come from the algorithm used. First of all, the "real curve" is continuous but the curve we work on is of course discrete. Phi is known on the points of the grid (it is calculated on the grid) and to build the curve, the algorithm used takes two neighbour points A, B of the grid (a size of a grid rectangle) and if the values of phi have an opposite signe then, as phi is assumed to be continuous, it means that there is a point on the segment between these two points where the value of phi is zero. But we know that the curve is defined as the set of points where phi equal zero. So the algorithm chooses a point on this segment with a linear rule between A and B to determine (knowing phi(A) and phi(B)) approximately where M, such as phi(M)=0, is). So we can notice that a point of the curve will be on a size of a grid rectangle. And the continuity of phi enables us to conclude either that two close points on a curve belong to a same rectangle (there can't be any "jumps" since it is continuous). The property used here is a maths results on continuous functions: the intermediate values theorem.

That's why, if two points are really close (distance less than the diagonal), we join them otherwise we don't.

Let s be the integer argument.

Note: We have seen that a single curve can be composed of different pieces.

Note: a piece is a close circuit in a curve. Example of a curve composed of several parts:

These following Matlab functions respectively display the curve number s for the first one, every curve which have a number N = 1 + s*k where N<total number of saved curves and k is an integer (s=1 prints all the saved curves) for the second one and the third one does exactely the same thing than 2 but this time the output is in avi format (a video) and it is saved in a file given as an argument. The first one and the second one have 1 argument: s and the third one takes 2 arguments: s and the name of the output file.

Source code of printseveral and explanations

function z = printseveral(s)

% N = PRINTCURVE(I)
%
% S: ...
%
% N: ...
%
% PRINT displays fronts saved in 'Curves'. ...

Loads all the files in matlab variables and initializes all the other variables	
  global x y n Fronts Time Phi

  Time = load('Time');
  
  Phi = load('Phi');
  
  Curves = load('Curves');
  CurveLengths = load('CurveLengths');
  x = load('X');
  y = load('Y');

  Delta_x = x(2) - x(1);
  Delta_y = y(2) - y(1);
  LimitDistance = sqrt( Delta_x^2 + Delta_y^2 );
  
  n = length(CurveLengths);

  Fronts = cell(n,1);

  count = 1;

Sort the curves and fill the array (more precisely a cell) Front with them	
  for i=1:n
    Fronts(i) = {Curves(count:count + CurveLengths(i) - 1, :)};
    count = count + CurveLengths(i);
  end

  figure;
  axis([x(1) x(end) y(1) y(end)]);

  hold on

Select the curve s	
   for i=1:s:n  

Initialize LastX & LastY.  LastX & LastY contain the first point of each new piece of the curve	   
    LastX = Fronts{i}(1, 1);
    LastY = Fronts{i}(1, 2);
 
Process the job with all the points except the "last one", the last one has to be linked with
 the first one of the last part of the curve (which can be in several pieces cf previously)
	for j=1:size(Fronts{i},1)-1

If distance<diagonal link the points
if ( sqrt( (Fronts{i}(j,1) - Fronts{i}(j+1,1))^2 + (Fronts{i}(j,2) - Fronts{i}(j+1,2))^2 )
 < LimitDistance)
plot([Fronts{i}(j, 1); Fronts{i}(j+1, 1)], [Fronts{i}(j, 2); Fronts{i}(j+1, 2)], '-b')

Else, close the current curve piece and store in LastX & LastY the first point of the new
 piece which is our second point
      else
	if ( sqrt( (Fronts{i}(j,1) - LastX)^2 + (Fronts{i}(j,2) - LastY)^2 ) < LimitDistance)
	  plot([Fronts{i}(j, 1); LastX], [Fronts{i}(j, 2); LastY], '-b')
	end
	LastX = Fronts{i}(j+1, 1);
	LastY = Fronts{i}(j+1, 2);
      end
    end

Processes the job with the last point which is linked to the first point of the last piece
    j = size(Fronts{i},1);
if ( sqrt( (Fronts{i}(j,1) - Fronts{i}(1,1))^2 + (Fronts{i}(j,2) - Fronts{i}(1,2))^2 )
 < LimitDistance)
      plot([Fronts{i}(j, 1); Fronts{i}(1, 1)], [Fronts{i}(j, 2); Fronts{i}(1, 2)], '-b')
    else
    if ( sqrt( (Fronts{i}(j,1) - LastX)^2 + (Fronts{i}(j,2) - LastY)^2 ) < LimitDistance)
	plot([Fronts{i}(j, 1); LastX], [Fronts{i}(j, 2); LastY], '-b')
      end
    end
  end
  set(gcf, 'Position', [100 100 600 600]);
  
  return;