6D-Reconstruction: added new default choices for data fidelity term and for TV regularization

data-fidelity: l_2 -> KL
TV: l_1 -> l_{1,2}

But the old ones remain available, even if they are not accessible to the user

Signed-off-by: Nicola Vigano <nicola.vigano@esrf.fr>

zUtil_Deformation/Gt6DBlobReconstructor.m

@@ -18,6 +18,9 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
 
         verbose = false;
 
+        detector_norm = 'KL'; % Possibilities are: {'KL'} | 'l2'
+        tv_norm = 'l12'; % Possibilities are: {'l12'} | 'l1'
+
         algo_ops_c_functions = true;
     end
 
@@ -317,15 +320,28 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
             dsES = gtMathsGradient(sES);
             self.statistics.add_timestamp(toc(c), 'cp_dual_update_tv', 'cp_dual_tv_gradient');
 
-            if (self.algo_ops_c_functions)
-                % Using the l1 update function because TV is the l1 of the
-                % gradient
-                q = gt6DUpdateDualL1_c(q, dsES, sigma, self.num_threads);
-            else
-                for n = 1:numel(q)
-                    q{n} = q{n} + sigma * dsES{n};
-                    q{n} = q{n} ./ max(1, abs(q{n}));
-                end
+            switch (self.tv_norm)
+                case 'l12'
+                    for n = 1:numel(q)
+                        q{n} = q{n} + sigma * dsES{n};
+                    end
+
+                    grad_l2 = sqrt(q{1} .^ 2 + q{2} .^ 2 + q{3} .^ 2);
+
+                    for n = 1:numel(q)
+                        q{n} = q{n} ./ max(1, grad_l2);
+                    end
+                case 'l1'
+                    if (self.algo_ops_c_functions)
+                        % Using the l1 update function because TV is the l1 of the
+                        % gradient
+                        q = gt6DUpdateDualL1_c(q, dsES, sigma, self.num_threads);
+                    else
+                        for n = 1:numel(q)
+                            q{n} = q{n} + sigma * dsES{n};
+                            q{n} = q{n} ./ max(1, abs(q{n}));
+                        end
+                    end
             end
 
             c = tic();
@@ -334,19 +350,30 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
         end
 
         function p = update_dual_detector(self, p, proj_bls, sigma, sigma_1)
-            if (self.algo_ops_c_functions)
-                p = gt6DUpdateDualDetector_c(p, self.blobs, proj_bls, sigma, sigma_1, self.num_threads);
-
-                % % Or equivalently
-                % proj_bls = internal_cell_sub_assign(proj_bls, bls);
-                % proj_bls = internal_cell_prod_assign(proj_bls, sigma1);
-                % p = internal_cell_sum_assign(p, proj_bls);
-                % p = internal_cell_prod_assign(p, sigma1_1);
-            else
-                for n = 1:numel(p)
-                    p{n} = p{n} + sigma{n} .* (proj_bls{n} - self.blobs{n});
-                    p{n} = p{n} .* sigma_1{n};
-                end
+            switch (self.detector_norm)
+                case 'KL'
+                    for n = 1:numel(p)
+                        temp_p = p{n} + sigma{n} .* proj_bls{n};
+                        temp_bls = sigma{n} .* self.blobs{n};
+                        temp_d = (temp_p - 1) .^ 2 + 4 .* temp_bls;
+                        temp_d = sqrt(abs(temp_d));
+                        p{n} = (1 + temp_p - temp_d) / 2;
+                    end
+                case 'l2'
+                    if (self.algo_ops_c_functions)
+                        p = gt6DUpdateDualDetector_c(p, self.blobs, proj_bls, sigma, sigma_1, self.num_threads);
+
+                        % % Or equivalently
+                        % proj_bls = internal_cell_sub_assign(proj_bls, bls);
+                        % proj_bls = internal_cell_prod_assign(proj_bls, sigma1);
+                        % p = internal_cell_sum_assign(p, proj_bls);
+                        % p = internal_cell_prod_assign(p, sigma1_1);
+                    else
+                        for n = 1:numel(p)
+                            p{n} = p{n} + sigma{n} .* (proj_bls{n} - self.blobs{n});
+                            p{n} = p{n} .* sigma_1{n};
+                        end
+                    end
             end
         end