%% test_zoom_3D.m
%
% E. Koenig, 08/04/2010
%
% Augmentation of information of a scale invariant 3D field by a factor N.
% The initial 3D field is created by a compound Poisson cascade (CPC)
% fractionally integrated. Its scale invariance properties are known as
%
% $$ \zeta(q) = \tau(q) + qH $$
%
% The variance and mean of this initial 3D field are modified by the
% parameters alpha and I_0 respectively. The initial field is of the form
%
% $$ I = I_0 + \alpha\mathcal{I}_{H}(Q-<Q>) $$
%
% where I is the initial field, Q is the CPC, < Q > is the mean of the
% cascade and the operator I_H stands for the pseudo-fractional 
% integration of parameter H.
%
% The augmentation of information by a factor N is performed as the
% succession of log2(N) elementary augmentations by a factor 2.
%
%
% This test code uses the functions:
%
% * <../../Functions/to_compile/CPC_3D/html/cpc_3D.html cpc_3D> to be compiled, in the folder CPC_3D
% * <../../Functions/html/zoom_x2_3D.html zoom_x2_3D>
%
%

%% Main code
%

%% Parameters
%
% parameters used to create the initial field
X = 31;     % we choose to create a cubic 3D field of size X^3
r_min = 1/X;    % finest resolution of the 3D field
dx = 1/X;   % sampling rate, this is the same in every direction
L = 1;  % highest resolution of the 3D field
law = 3;  % distribution law used to impose the desired scale 
          % invariance properties of the 3D field
param = [0.1 0 0.5]; % law parameters
pattern = 2;  % pattern of the field
angle = 0;  % use of a non variable angle
seed_in = 0;  % initial seed used for the random drafts

H = 0.7;    % parameter for the pseudo-fractional integration used to smooth
            % the field

alpha_ref = 8;  % parameter used to modify the variance of the field
I_0_ref = 23;   % parameter used to modify the mean of the field

N = 2;  % number of successive elementary augmentations by a factor 2



%% Synthesis of scale invariant 3D field
%
% creation of the compound Poisson cascade (CPC) with prescribed scale invariance
% properties
[Q seed_out] = cpc_3D(L, L, L, r_min, dx, law, param, pattern, angle, seed_in);
seed_in = seed_out;
Q = reshape(Q, X+1, X+1, X+1);
% modification of the CPC and pseudo-fractional integration
Q = alpha_ref.*(Q-mean(Q(:)));
I_0 = I_0_ref + integfrac3D(Q, H);

clear Q


%% Augmentation procedure

tic;

for n = 1:N
    fprintf(1,'Augmentation of the initial field by a factor %d\n', 2^n);
    eval(['[I_',int2str(n),', seed_out] = zoom_x2_3D(I_',int2str(n-1),', law, param, H, 0, 0, seed_in);']);
    seed_in = seed_out;
end

toc



%% Visualisation using Matlab

fprintf(1,'Visualisation using Matlab\n');

for n = 0:N
    eval(['[x y z] = meshgrid(1:size(I_',int2str(n),',1), 1:size(I_',int2str(n),',2), 1:size(I_',int2str(n),',3));']);
    figure
    eval(['p=patch(isosurface(x, y, z, I_',int2str(n),', mean(I_0(:))/2^(3*n)));']);
    eval(['isonormals(x, y, z, I_',int2str(n),', p);']);
    set(p,'FaceColor','red','EdgeColor','none');
    daspect([1,1,1]);
    view(3);
    axis tight
    camlight
    lighting phong
    eval(['print -dpng images/I_', int2str(n),'_Matlab']);
end



%% Script for a visualisation using Marching Cubes
%
% Write the correct .dat file to be used for a visualization using Marching
% Cubes.
% 
% To use our visualization method :
%
%   cd ../Visualization/
%   ./champ_Mac ../Test/images/I_1.dat 10
%

fprintf(1,'Script for a visualisation using Marching Cubes\n');

for n=0:N
    eval(['I_write = (I_',int2str(n),'-min(I_',int2str(n),'(:)))/(max(I_',int2str(n),'(:))-min(I_',int2str(n),'(:)));']);
    eval(['f = fopen(''images/I_',int2str(n),'.dat'', ''wt'');']);
    eval(['S = size(I_',int2str(n),');']);
    fprintf(f, '%d %d %d\n', S(1), S(2), S(3));
    for x=1:S(1)
      for y=1:S(2)
          for z=1:S(3)
              fprintf(f, '%f\n', I_write(x,y,z));
          end
      end
    end
    fclose(f);
end


