%% zoom_x2_3D.m
%
% E. Koenig, 08/04/2010
%
% Augmentation of information of a 3D field by a factor 2.
%
%
%% References
%
% |[1] Koenig, E., Chainais, P.: Virtual resolution enhancement of scale 
% invariant textured images using stochastic processes, Proceedings of ICIP 
% 2009 (2009)|
%
% |[2] Chainais, P.: Infinitely divisible cascades to model the statistics 
% of natural images, IEEE Trans. on P.A.M.I. 29(12), 2105-2119 (2007)|
%

%% Description
%
% The size of the initial field is multiplied by two while new small scale
% information is added below the initial voxel resolution.
% The initial 3D field is supposed to be scale invariant with known scale
% invariance parameters.
% Furthermore, it can be modelized by a 3D compound Poisson cascade (CPC)
% [2]. 
%
% This field is seen as an integrated 3D CPC modified in variance by alpha 
% and in mean by I_0:
%
% $$ I = I_0 + \alpha\mathcal{I}_{H}(Q-<Q>) $$
%
%
% The high resolution field must have the same scale invariance properties
% as the initial field but on a larger range of scales.
% The new information added must be scale invariant. 
% It is created by the stochastic model of CPC [2]. 
%
%
% The augmentation method is divided into 5 parts:
%
% * *interpolation*: interpolates the initial field to the new resolution.
% This step adds no information to the initial field. It is realized using
% splines.
% * *fractional derivation*: apply a fractional derivation to access the
% underlying CPC (Q) of the initial field.
% * *augmentation of information*: add new details (created by CPC) with
% the same scale invariant properties as the initial field.
% * *fractional integration*: apply a fractional integration to get the new
% augmentaed field with the right scale invariant properties.
% * *conservation*: ensures that a zoom out to the initial resolution gives
% the initial field back.
%
%
% For more information on this method see [1].
%
%
% Needs the functions:
%
% * <interp_x2_3D.html interp_x2_3D>
% * <derivfrac3D.html derivfrac3D>
% * <givexpIDC2D.html givexpIDC2D>
% * <../to_compile/CPC_3D/html/cpc_3D.html cpc_3D> to be compiled, in the folder CPC_3D
% * <integfrac3D.html integfrac3D>
% * <zoom_out_x2_3D.html zoom_out_x2_3D>
% 

%% Use
%
% [I_aug, seed_out] = zoom_x2_3D(I, law, param, H, alpha_ref, I_0_ref, seed_in);
%

%% Inputs
%
% * *I*:  input 3D field
% * *law*: distribution law used to create the details with good
% properties: 
%
%          - 1: log-normal
%          - 2: exponential
%          - 3: log-Poisson
%
% * *param*: parameters of the law. The law and its parameters are 
% determined by a previous analysis of the scale invariance properties of 
% the initial field.
% * *H*: parameter of the fractional derivation and integration.
% * *alpha_ref*: if known, this is the value of the modified variance of
% the initial field. 0 if not known.
% * *I_0_ref*: if known, this is the value of the modified mean of the 
% initial field. 0 if not known.
% * *seed_in*: value of the initial seed used to create the stochastic
% details.
%

%% Outputs
%
% * *I_aug*: augmented 3D field. This field has the same scale invariance
% properties as the initial field but on a larger range of scales, and a
% zoom out to the initial resolution yields the initial field back.
% * *seed_out*: last seed used during the process. It can be used to
% iterate the augmentation process whithout having redundances when
% creating small scale details.
%

%% Example
%
% For a complete example of use, see the test program <../../Test/html/test_zoom_3D.html test_zoom_3D>.
%
% 


function [I_aug, seed_out] = zoom_x2_3D(I, law, param, H, alpha_ref, I_0_ref, seed_in)

% size of the input field
S = size(I);
% zoom factor, power of 2
zoomfactor = 2;

S_aug = S.*zoomfactor;


% Conservative interpolation -> I_interp
% based on the works of Thevenaz
I_int = interp_x2_3D(I, '*splines', 0);
I_interp = max(0,I_int);    %positivity constraint on I_interp

clear I_int

% Pseudo-fractionnal integration -> DH_I
% x||k||^H in the Fourier domain
DH_I = derivfrac3D(I_interp,H);

%  Estimation of alpha and I_0
q = [-2:0.5:0 1/6:1/6:5];
tau_theo = givexpIDC2D(q, law, param);
tau2 = tau_theo(17);
r0 = 1/max(S);
if(alpha_ref~=0)
    alpha = alpha_ref;
else
    alpha = sqrt(var(DH_I(:))/((r0^tau2)-1));
end
if(I_0_ref~=0)
    I_0 = I_0_ref;
else
    I_0 = mean(I_interp(:));
end

clear alpha_ref I_0_ref I_interp tau2 tau_theo q

% Computation of small scales using CPC -> Qr_r2
% r = original resolution
r = r0;
% r2 = new resolution
r2 = r/zoomfactor;
% maximum scale
xmax = max(S);
% use of cpc_3D
[Q seed_out] = cpc_3D(xmax, xmax, xmax, r2/r, r2/r, law, param, 2, 0, seed_in);
Qr_r2 = reshape(Q, max(S_aug)+1, max(S_aug)+1, max(S_aug)+1);
Qr_r2 = Qr_r2(1:S_aug(1), 1:S_aug(2), 1:S_aug(3));

clear Q r0 r2 xmax r S_aug S zoomfactor law param


% Combination of high and small scales -> DH_J
DH_J = DH_I.*Qr_r2 + alpha.*1.*(Qr_r2-mean(Qr_r2(:)));

clear DH_I Qr_r2 alpha

% Pseud-fractionnal integration -> J
% /||k||^H in the Fourier domain
J = integfrac3D(DH_J,H);

clear DH_J H

% Add the mean 
J = J+I_0;
clear I_0
% if J is negative, there is a problem
% a natural field can't be negative (accumulation of photons)
% I_aug null
count_neg = 0;
if(min(J(:))<0)
    neg_int = J(J(:)<=0);
    find(J(:)<=0);
    count_neg = length(neg_int);
    % solution 1 : saturation to 0.0
    J(J(:)<=0) = 0;
end
clear neg_int


% Conservation -> I_aug
% backward zoom
J_int = zoom_out_x2_3D(J);
J_conv_dev = J_int;

clear J_int

% conservation map: ratio between the new and the input fields
J_norm = zeros(size(I));
listok = find(I(:)>0 & J_conv_dev(:)>0);
J_norm(listok) = I(listok) ./ J_conv_dev(listok);

clear J_conv_dev I listok

% interpolation of the map to the new resolution
J_int = interp_x2_3D(J_norm,'*nearest');
J_norm2 = J_int;

clear J_norm J_int pas i

% conservation of J and creation of the augmented field
I_aug = J .* J_norm2;


clear J_norm2 J

