attended to reviewers remarks on computing attenuation factors

src/xSTIR/cSTIR/cstir.cpp

@@ -650,21 +650,21 @@ void* cSTIR_createPETAttenuationModel(const void* ptr_img, const void* ptr_am)
 }
 
 extern "C"
-void* cSTIR_computeACF(const void* ptr_sino, const void* ptr_mumap, const void* ptr_asm,
-    void* ptr_att, void* ptr_acf, void* ptr_iacf)
+void* cSTIR_computeACF(const void* ptr_sino, const void* ptr_mumap, const void* ptr_am,
+    const void* ptr_att, void* ptr_af, void* ptr_acf)
 {
 	try {
 		SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_sino, ptr_sino);
 		SPTR_FROM_HANDLE(STIRImageData, sptr_mumap, ptr_mumap);
+		SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_af, ptr_af);
 		SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_acf, ptr_acf);
-		SPTR_FROM_HANDLE(STIRAcquisitionData, sptr_iacf, ptr_iacf);
-		SPTR_FROM_HANDLE(AcqMod3DF, sptr_asm, ptr_asm);
+		SPTR_FROM_HANDLE(AcqMod3DF, sptr_am, ptr_am);
 		SPTR_FROM_HANDLE(PETAttenuationModel, sptr_att, ptr_att);
 		PETAttenuationModel::compute_ac_factors(sptr_sino, sptr_mumap,
-			sptr_asm, sptr_att, sptr_acf, sptr_iacf);
+			sptr_am, sptr_att, sptr_af, sptr_acf);
+		HANDLE_FROM_SPTR(STIRAcquisitionData, sptr_af, ptr_af);
 		HANDLE_FROM_SPTR(STIRAcquisitionData, sptr_acf, ptr_acf);
-		HANDLE_FROM_SPTR(STIRAcquisitionData, sptr_iacf, ptr_iacf);
-		HANDLE_FROM_SPTR(PETAttenuationModel, sptr_att, ptr_att);
+		//HANDLE_FROM_SPTR(PETAttenuationModel, sptr_att, ptr_att);
 		return (void*) new DataHandle;
 	}
 	CATCH;

src/xSTIR/cSTIR/include/sirf/STIR/cstir.h

@@ -79,7 +79,7 @@ extern "C" {
 	void* cSTIR_createPETAttenuationModel
 		(const void* ptr_img, const void* ptr_am);
 	void* cSTIR_computeACF
-		(const void* ptr_sino, const void* ptr_mumap, const void* ptr_am, void* ptr_att, void* ptr_acf, void* ptr_iacf);
+		(const void* ptr_sino, const void* ptr_mumap, const void* ptr_am, const void* ptr_att, void* ptr_acf, void* ptr_iacf);
 	void* cSTIR_chainPETAcquisitionSensitivityModels
 		(const void* ptr_first, const void* ptr_second);
 	void* cSTIR_setupAcquisitionSensitivityModel(void* ptr_sm, void* ptr_ad);

src/xSTIR/cSTIR/include/sirf/STIR/stir_x.h

@@ -552,28 +552,31 @@ The actual algorithm is described in
 		virtual void unnormalise(STIRAcquisitionData& ad) const;
 		// divide by bin efficiencies (here attenuation factors), i.e. correct data in \a ad for attenuatio
 		virtual void normalise(STIRAcquisitionData& ad) const;
+		/*! Convenience function computing attenuation factor using unnormalise
+		    and its inverse (attenuation correction factor) using STIRImageData::inv
+		*/
 		static void compute_ac_factors(
+			// input arguments
 			std::shared_ptr<STIRAcquisitionData> acq_templ_sptr,
 			std::shared_ptr<STIRImageData> mu_map_sptr,
 			std::shared_ptr<PETAcquisitionModel>& am_sptr,
 			std::shared_ptr<PETAttenuationModel>& asm_sptr,
-			std::shared_ptr<STIRAcquisitionData>& acf_sptr,
-			std::shared_ptr<STIRAcquisitionData>& iacf_sptr)
+			// output arguments
+			std::shared_ptr<STIRAcquisitionData>& af_sptr,
+			std::shared_ptr<STIRAcquisitionData>& acf_sptr)
 		{
 			//PETAcquisitionModelUsingRayTracingMatrix acq_mod;
 			PETAcquisitionModel& acq_mod = *am_sptr;
 			acq_mod.set_up(acq_templ_sptr, mu_map_sptr);
-			asm_sptr = std::shared_ptr<PETAttenuationModel>(new PETAttenuationModel(*mu_map_sptr, acq_mod));
-//			std::shared_ptr<PETAttenuationModel>
-//				asm_sptr(new PETAttenuationModel(*mu_map_sptr, acq_mod));
 			PETAttenuationModel& acq_sens_mod = *asm_sptr;
 			acq_sens_mod.set_up(acq_templ_sptr->get_exam_info_sptr(),
 				acq_templ_sptr->get_proj_data_info_sptr()->create_shared_clone());
-			acf_sptr = acq_templ_sptr->clone();
-			acf_sptr->fill(1.0);
-			iacf_sptr = acf_sptr->clone();
-			acq_sens_mod.unnormalise(*acf_sptr);
-			acq_sens_mod.normalise(*iacf_sptr);
+			af_sptr = acq_templ_sptr->clone();
+			af_sptr->fill(1.0);
+			acf_sptr = af_sptr->clone();
+			acq_sens_mod.unnormalise(*af_sptr);
+			//acq_sens_mod.normalise(*acf_sptr);
+			acf_sptr->inv(0, *af_sptr);
 		}
 
 	protected:

src/xSTIR/cSTIR/tests/test7.cpp

@@ -88,12 +88,12 @@ int main()
 		std::cout << "===== sinograms norm: " << sinograms_sptr->norm() << '\n';
 		std::cout << "===== randoms norm: " << randoms_sptr->norm() << '\n';
 
-		std::shared_ptr<PETAttenuationModel> att_sptr;
-		std::shared_ptr<STIRAcquisitionData> acf_sptr; // attenuation correction factor
-		std::shared_ptr<STIRAcquisitionData> iacf_sptr; // the inverse of the above
-		PETAttenuationModel::compute_ac_factors(sinograms_sptr, mu_map_sptr, am_sptr, att_sptr, acf_sptr, iacf_sptr);
-		std::cout << acf_sptr->norm() << '\n';
-		std::cout << iacf_sptr->norm() << '\n';
+		std::shared_ptr<PETAttenuationModel> att_sptr = std::shared_ptr<PETAttenuationModel>(new PETAttenuationModel(mu_map, am));
+		std::shared_ptr<STIRAcquisitionData> af_sptr; // attenuation factor
+		std::shared_ptr<STIRAcquisitionData> acf_sptr; // the inverse of the above
+		PETAttenuationModel::compute_ac_factors(sinograms_sptr, mu_map_sptr, am_sptr, att_sptr, af_sptr, acf_sptr);
+		std::cout << "===== norm of the attenuation factor: " << af_sptr->norm() << '\n';
+		std::cout << "===== norm of the attenuation correction factor: " << acf_sptr->norm() << '\n';
 
 		CREATE_OBJ(PETAcquisitionSensitivityModel, acq_sm, acq_sm_sptr, f_norm);
 
@@ -102,9 +102,9 @@ int main()
 		se.set_attenuation_image_sptr(mu_map_sptr);
 		se.set_background_sptr(randoms_sptr);
 		se.set_asm(acq_sm_sptr);
-		se.set_attenuation_correction_factors_sptr(iacf_sptr);
+		se.set_attenuation_correction_factors_sptr(acf_sptr);
 		se.set_num_iterations(4);
-		std::cout << "number of scatter iterations that will be used: " << se.get_num_iterations() << '\n';
+		std::cout << "===== number of scatter iterations that will be used: " << se.get_num_iterations() << '\n';
 		se.set_OSEM_num_subsets(7);
 		se.set_output_prefix("scatter");
 		se.set_up();
@@ -114,7 +114,7 @@ int main()
 		std::cout << "===== scatter estimate norm: " << scatter_sptr->norm() << '\n';
 
 		acq_sm.set_up(*acf_sptr);
-		std::shared_ptr<STIRAcquisitionData> mf_sptr = acf_sptr->clone();
+		std::shared_ptr<STIRAcquisitionData> mf_sptr = af_sptr->clone();
 		acq_sm.unnormalise(*mf_sptr);
 		std::cout << "===== multfactors norm: " << mf_sptr->norm() << '\n';
 

src/xSTIR/pSTIR/STIR.py

@@ -1740,12 +1740,12 @@ def compute_attenuation_factors(sinograms, mu_map):
         '''Creates attenuation model and returns it as an AcquisitionSensitivityModel object
            together with the attenuation factor and its inverse as AcquisitionData objects
         '''
-        asm = AcquisitionModelUsingRayTracingMatrix()
-        attn = AcquisitionSensitivityModel(mu_map, asm)
+        am = AcquisitionModelUsingRayTracingMatrix()
+        attn = AcquisitionSensitivityModel(mu_map, am)
+        af = AcquisitionData(sinograms)
         acf = AcquisitionData(sinograms)
-        iacf = AcquisitionData(sinograms)
-        try_calling(pystir.cSTIR_computeACF(sinograms.handle, mu_map.handle, asm.handle, attn.handle, acf.handle, iacf.handle))
-        return attn, acf, iacf
+        try_calling(pystir.cSTIR_computeACF(sinograms.handle, mu_map.handle, am.handle, attn.handle, af.handle, acf.handle))
+        return af, acf
 
     def __del__(self):
         """del."""