Fix superior dose limit in target

src/extra/opt4D_planfile_template.pln

@@ -30,11 +30,12 @@ OBJ 1  dose_file_name missing_target_dose.dat
 OBJ 1  min_dose 20000
 OBJ 1  weight 30
 
+# Create upper limit at target dose+3 Gy (60+3 Gy)
 OBJ 2  type voxeldose
 OBJ 2  VOI 1
 OBJ 2  dose_file_name missing_target_dose.dat
 OBJ 2  max_dose 20000
-OBJ 2  scale_factor 1.025
+OBJ 2  sum_constant 3
 OBJ 2  weight 10
 
 VOI 2 name falloff

src/program_options.cpp

@@ -297,8 +297,8 @@ void correct_energy_range (const RShifter_steps_t& rshifter_steps,
                 if (no_errors) {
                     std::copy(temp.begin(), temp.end(), energies.begin() + accu_spots);
                     // Create range shifter to make it deliverable
-                    std::cerr << "New range shifter for beam " << ibeam << " = ";
-                    std::cerr << rshifter << " water cm." << std::endl;
+                    std::cout << "New range shifter for beam " << ibeam << " = ";
+                    std::cout << rshifter << " water cm." << std::endl;
                     if (rs.at(ibeam).exists) {
                         rs.at(ibeam).set_wepl(rshifter);
                     } else {

src/utils/opt4D_manager.cpp

@@ -351,10 +351,8 @@ void Opt4D_manager::set_write_reference_dose(const std::vector<float>& dose)
     std::vector<float> ave_ref_batches(nbatches, 0);
     float ave_ref = 0;
     for(size_t i = 0; i<n_voxels; i++) {
-        if (min_ref > dose.at(i))
-            min_ref = dose.at(i);
-        if (max_ref < dose.at(i))
-            max_ref = dose.at(i);
+        min_ref = std::min(min_ref, dose.at(i));
+        max_ref = std::max(max_ref, dose.at(i));
         ave_ref_batches.at(int(i%nbatches)) += dose.at(i)/n_spots;
     }
     for (uint i = 0; i < nbatches; ++i) {

src/weight_adaptation.cpp

@@ -271,6 +271,11 @@ void adapt_weights (const std::vector<std::string>& dij_files,         // Dij fi
             dose_by_subset[i] += d;
             baseline_dose[i] -= d;
             if (target_mask.data.at(i) > 0.5) {
+                // the result could be negative!! But that would be fine. Negative voxels are set to
+                // zero by Opt4D. The reason not to do it here is that I have allowed Opt4D to add a
+                // constant value to the voxel-dose to control the maximum dose within the target.
+                // This operation needs to know if the voxel is -3 or -10 or whatever value. If
+                // after summing, the value is still negative, it will be set to zero.
                 target_dose_missing[idij] = dose_prescription - baseline_dose[i];
             }
             idij++;