%% cpc_3D.m
%
% E. Koenig, 08/04/2010
%
% Synthesis of a scale invariant 3D field using compound Poisson cascades
% [1].
%
% The field is characterized by scale invariant exponents 
%
% $$ \tau(q) $$
%
% prescribed by the distribution law of the values of random i.i.d. positive 
% variables.
% This scale invariance is obtained in the range of scales between L and l
% where L is the largest scale of the field and l is the smallest one.
%
%
%% Reference
%
% |[1] Chainais, P.: Infinitely divisible cascades to model the statistics 
% of natural images, IEEE Trans. on P.A.M.I. 29(12), 2105-2119 (2007)|
%
%
%
%% Use
%
% [Q seed_out] = cpc_3D(L_x, L_y, L_z, r_min, dx, law, param, pattern, angle, seed_in);
%
%

%% Inputs
%
% * *L_x, L_y, L_z*: largest scales of the cascade in each dimension
% * *r_min*: smallest scale of the cascade 
% * *dx*: sampling rate, this is the same in every direction
% * *law*: distribution law used to impose the desired scale invariance 
% properties of the 3D field
%
%          - 1: log-normal
%          - 2: exponential
%          - 3: log-Poisson
%
% * *param*: law parameters
% * *pattern*: pattern of the field
%
%          - 1: cubes
%          - 2: balls
%          - 3: shaper square cosine
%          - 4: shaper exponential
%
% * *angle*: 1 - variable angle, 0 - fix angle
% * *seed_in*: value of the initial seed used for the various random
% variables to draft
%
%

%% Outputs
%
% * *Qr*: 3D resulting field. The output is in 1 dimension, it needs a
% reshape using the |reshape| function of Matlab
% * *seed_out*: last seed used by the random variable generator
%
%

%% Example
%
% Creation of a 64x64x64 scale invariant 3D field
%

% Parameters
X = 63;
L_x = 4;
L_y = 4;
L_z = 4;
r_min = 4/16;
dx = 4/X;
law = 2;
param = [0.7 0 1];
pattern = 2;
angle = 0;
seed_in = 0;

% Computation of the 3D CPC
[Q seed_out] = cpc_3D(L_x, L_y, L_z, r_min, dx, law, param, pattern, angle, seed_in);
Q = reshape(Q, X+1, X+1, X+1);

% Visualisation using Matlab
[x y z] = meshgrid(1:size(Q,1), 1:size(Q,2), 1:size(Q,3));

figure
p=patch(isosurface(x, y, z, Q, mean(Q(:))));
isonormals(x, y, z, Q, p);
set(p,'FaceColor','red','EdgeColor','none');
daspect([1,1,1]);
view(3);
axis tight
camlight
lighting phong

%%