6D-reconstruction: refactored use of lambdas

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

5_reconstruction/gtReconstruct6DLaunchAlgorithm.m

@@ -25,6 +25,9 @@ function algo = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, parameters, va
     algo.tv_norm = rec_opts.tv_norm;
     algo.tv_strategy = rec_opts.tv_strategy;
 
+    algo.lambda_l1 = rec_opts.lambda_l1;
+    algo.lambda_tv = rec_opts.lambda_tv;
+
     geom_memory = sum(arrayfun(@(x)x.get_memory_consumption_geometry(), sampler));
     grain_memory = sum(arrayfun(@(x)x.get_memory_consumption_graindata(), sampler));
 
@@ -39,5 +42,5 @@ function algo = gtReconstruct6DLaunchAlgorithm(sampler, rec_opts, parameters, va
         'Sinograms %d MB, Geometry %d MB, Grain data %d MB)\n'], ...
         round(mem_cons / (2^20)))
 
-    algo.cp(upper(rec_opts.algorithm), rec_opts.num_iter, 'lambda_l1', rec_opts.lambda_l1, 'lambda_tv', rec_opts.lambda_tv);
+    algo.cp(upper(rec_opts.algorithm), rec_opts.num_iter);
 end

zUtil_Deformation/Gt6DBlobReconstructor.m

@@ -19,6 +19,9 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
         tv_norm = Gt6DBlobReconstructor.possible_tv_norms{1};
         tv_strategy = Gt6DBlobReconstructor.possible_tv_strategies{1};
 
+        lambda_l1 = 1e-2;
+        lambda_tv = 1;
+
         orientation_groups = [];
 
         algo_ops_c_functions = true;
@@ -163,21 +166,19 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
             self.cp('6DLS', numIters);
         end
 
-        function cp_l1(self, numIters, lambda)
-            self.cp('6DL1', numIters, 'lambda_l1', lambda)
+        function cp_l1(self, numIters)
+            self.cp('6DL1', numIters)
         end
 
-        function cp_tv(self, numIters, lambda)
-            self.cp('6DTV', numIters, 'lambda_tv', lambda);
+        function cp_tv(self, numIters)
+            self.cp('6DTV', numIters);
         end
 
-        function cp_tvl1(self, numIters, lambda)
-            self.cp('6DTVL1', numIters, lambda);
+        function cp_tvl1(self, numIters)
+            self.cp('6DTVL1', numIters);
         end
 
-        function cp(self, algo, numIters, varargin)
-            conf = struct('lambda_l1', [], 'lambda_tv', []);
-            conf = parse_pv_pairs(conf, varargin);
+        function cp(self, algo, numIters)
 
             self.statistics.clear();
 
@@ -198,10 +199,10 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
             fprintf('Reconstruction using algorithm: %s\n', upper(algo));
             fprintf(' - Detector data-fidelity term: %s\n', self.detector_norm);
             if (do_l1_update)
-                fprintf(' - l1-term lambda: %g\n', conf.lambda_l1);
+                fprintf(' - l1-term lambda: %g\n', self.lambda_l1);
             end
             if (do_tv_update)
-                fprintf(' - TV-term lambda: %g\n', conf.lambda_tv);
+                fprintf(' - TV-term lambda: %g\n', self.lambda_tv);
                 fprintf(' - TV norm: %s\n', self.tv_norm);
                 fprintf(' - TV strategy: %s\n', self.tv_strategy);
             end
@@ -233,7 +234,7 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
                 % Computing update dual TV
                 if (do_tv_update)
                     self.statistics.tic('cp_dual_update_tv');
-                    [q_tv, mdiv_tv] = self.update_dual_TV(q_tv, nextEnhancedSolution, conf.lambda_tv);
+                    [q_tv, mdiv_tv] = self.update_dual_TV(q_tv, nextEnhancedSolution);
                     self.statistics.toc('cp_dual_update_tv');
                 end
 
@@ -251,13 +252,13 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
                         compute_update_primal = @(ii){};
                     case '6DL1'
                         compute_update_primal = @(ii)gtCxxMathsCellTimes(...
-                            q_l1(ii), conf.lambda_l1, 'threads', self.num_threads);
+                            q_l1(ii), self.lambda_l1, 'threads', self.num_threads);
                     case '6DTV'
                         compute_update_primal = @(ii)get_divergence(self, mdiv_tv, ii);
                     case '6DTVL1'
                         compute_update_primal = @(ii)gtCxxMathsCellPlus(...
                             get_divergence(self, mdiv_tv, ii), q_l1(ii), ...
-                            'scale', conf.lambda_l1, 'threads', self.num_threads);
+                            'scale', self.lambda_l1, 'threads', self.num_threads);
                 end
 
                 % Computing update primal
@@ -273,7 +274,7 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
 
                 if (self.verbose && (mod(ii, sample_rate) == 0))
                     self.normResiduals(ii / sample_rate) ...
-                        = self.compute_functional_value(p, algo, conf.lambda_l1);
+                        = self.compute_functional_value(p, algo);
                 end
                 drawnow();
             end
@@ -348,10 +349,7 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
     %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
     %%% Primal and Duals updating functions
     methods (Access = public)
-        function [q, mdivq] = update_dual_TV(self, q, new_enh_sol, lambda)
-            if (~exist('lambda', 'var') || isempty(lambda))
-                lambda = 1;
-            end
+        function [q, mdivq] = update_dual_TV(self, q, new_enh_sol)
 
             switch (self.tv_strategy)
                 case 'volume'
@@ -447,7 +445,7 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
 
             c = tic();
             for ii_g = 1:num_groups
-                mdivq{ii_g} = - lambda * gtMathsDivergence(q(ii_g, :));
+                mdivq{ii_g} = - self.lambda_tv * gtMathsDivergence(q(ii_g, :));
             end
             divergence_time = toc(c);
 
@@ -663,15 +661,15 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
                     end
                 case '6DL1'
                     for n = 1:num_geoms
-                        tau{n} = cast(- 1 ./ (self.bwd_weights{n} + 1), self.data_type);
+                        tau{n} = cast(- 1 ./ (self.bwd_weights{n} + 1 * self.lambda_l1), self.data_type);
                     end
                 case '6DTV'
                     for n = 1:num_geoms
-                        tau{n} = cast(- 1 ./ (self.bwd_weights{n} + 6), self.data_type);
+                        tau{n} = cast(- 1 ./ (self.bwd_weights{n} + 6 * self.lambda_tv), self.data_type);
                     end
                 case '6DTVL1'
                     for n = 1:num_geoms
-                        tau{n} = cast(- 1 ./ (self.bwd_weights{n} + 6 + 1), self.data_type);
+                        tau{n} = cast(- 1 ./ (self.bwd_weights{n} + 6 * self.lambda_tv + 1 * self.lambda_l1), self.data_type);
                     end
                 otherwise
             end
@@ -733,7 +731,7 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
             end
         end
 
-        function [value, summed_vols] = compute_functional_value(self, p, algo, lambda)
+        function [value, summed_vols] = compute_functional_value(self, p, algo)
             value = gtMathsDotProduct(p, self.blobs);
             num_geoms = self.get_number_geometries();
 
@@ -756,11 +754,11 @@ classdef Gt6DBlobReconstructor < Gt6DVolumeToBlobProjector
             end
             if (any(strcmpi(algo, {'6DTV', '6DTVL1'})))
                 sES = gtMathsGradient(summed_vols / num_geoms);
-                value = value + gtMathsNorm_l1(sES);
+                value = value + self.lambda_tv * gtMathsNorm_l1(sES);
             end
 
             if (any(strcmpi(algo, {'6DL1', '6DTVL1'})))
-                value = value + lambda * gtMathsNorm_l1(self.currentSolution);
+                value = value + self.lambda_l1 * gtMathsNorm_l1(self.currentSolution);
             end
         end