From 20de3a063ff5a7e603c15f985fd5708a8d374bbd Mon Sep 17 00:00:00 2001 From: TommasoFerri97 <79870251+TommasoFerri97@users.noreply.github.com> Date: Mon, 31 Aug 2026 15:17:32 +0200 Subject: [PATCH 1/3] Add channel edge handling to PinholeSPECTUB --- .../recon_buildblock/PinholeSPECTUB_Tools.h | 2 + src/recon_buildblock/PinholeSPECTUB_Tools.cxx | 17 +- .../PinholeSPECTUB_Weight3d.cxx | 362 ++++++++++++++++-- 3 files changed, 347 insertions(+), 34 deletions(-) diff --git a/src/include/stir/recon_buildblock/PinholeSPECTUB_Tools.h b/src/include/stir/recon_buildblock/PinholeSPECTUB_Tools.h index c552869f45..9f8b71081c 100644 --- a/src/include/stir/recon_buildblock/PinholeSPECTUB_Tools.h +++ b/src/include/stir/recon_buildblock/PinholeSPECTUB_Tools.h @@ -89,6 +89,7 @@ typedef struct bool do_round; // true: round shape || false: rectangular shape float dxcm; // horizontal size of the hole (cm): horizontal axis, diameter float dzcm; // vertical size of the hole (cm): vertical axis, diameter + float dycm; // channel edge width (cm) } hole_type; @@ -97,6 +98,7 @@ typedef struct typedef struct { std::string model; // cylindrical (cyl) or polygonal prism (pol) + std::string type; // knife-edge (knife) or channel-edge (channel) float rad; // radius of cylinder containig holes (cyl) or apothem (pol) float L; // collimator thickness diff --git a/src/recon_buildblock/PinholeSPECTUB_Tools.cxx b/src/recon_buildblock/PinholeSPECTUB_Tools.cxx index 103f3afc63..1570041699 100644 --- a/src/recon_buildblock/PinholeSPECTUB_Tools.cxx +++ b/src/recon_buildblock/PinholeSPECTUB_Tools.cxx @@ -442,6 +442,10 @@ read_coll_params_mph(wmh_mph_type& wmh) wmh.collim.model = wm_SPECT_read_value_1d(&stream1, DELIMITER); + // collimator type: knife or channel + token = wm_SPECT_read_value_1d(&stream1, DELIMITER); + wmh.collim.type = token; + token = wm_SPECT_read_value_1d(&stream1, DELIMITER); wmh.collim.rad = std::stof(token); @@ -481,7 +485,8 @@ read_coll_params_mph(wmh_mph_type& wmh) stream1.close(); info_stream << "Collimator parameters" << endl; - info_stream << "\nCollimator model: " << wmh.collim.model << endl; + info_stream << "\nCollimator type: " << wmh.collim.type << endl; + info_stream << "Collimator model: " << wmh.collim.model << endl; info_stream << "Collimator rad: " << wmh.collim.rad << endl; info_stream << "Number of holes: " << wmh.collim.Nht << endl; @@ -614,6 +619,16 @@ wm_SPECT_read_hvalues_mph(ifstream* stream1, char DELIMITER, int* nh, bool do_cy max_hsxcm = h.dxcm; pos1 = pos3; + //... channel-edge width y cm ....................... + + pos2 = line.find_first_not_of(" \t\f\v\n\r", pos1 + 1); + pos3 = line.find_first_of(" \t\f\v\n\r", pos2); + if (pos2 == string::npos) + error_wmtools_SPECT_mph(333, *nh, "dycm"); + token = line.substr(pos2, pos3 - pos2); + h.dycm = std::stof(token); + pos1 = pos3; + //... dimension z cm ....................... pos2 = line.find_first_not_of(" \t\f\v\n\r", pos1 + 1); diff --git a/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx b/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx index e3e8d590a5..b05655c023 100644 --- a/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx +++ b/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx @@ -474,59 +474,355 @@ check_zang_par(const voxel_type* v, const hole_type* h) void voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const wmh_mph_type& wmh) { + if (wmh.collim.type == "knife" || h->dycm == 0.) //if the channel width is set to zero the collimator is a knife-edge pinhole + { - //...vector voxel-hole, angles and distances............................... + //...vector voxel-hole, angles and distances............................... - float ux1 = h->x1 - v->x1; - float uy1 = h->y1 - v->y1; - float uz1 = h->z1 - v->z; + float ux1 = h->x1 - v->x1; + float uy1 = h->y1 - v->y1; + float uz1 = h->z1 - v->z; - if (uy1 <= EPSILON) - error_weight3d(88, ""); + if (uy1 <= EPSILON)error_weight3d(88, ""); - //...vector voxel-hole and distances............................... + //...vector voxel-hole and distances............................... - float dxyz_2 = ux1 * ux1 + uy1 * uy1 + uz1 * uz1; - float dvh_l = sqrtf(dxyz_2); + float dxyz_2 = ux1 * ux1 + uy1 * uy1 + uz1 * uz1; + float dvh_l = sqrtf(dxyz_2); - ux1 /= dvh_l; - uy1 /= dvh_l; - uz1 /= dvh_l; + ux1 /= dvh_l; + uy1 /= dvh_l; + uz1 /= dvh_l; - //...distance over voxel-hole line from voxel and hole to detection plane.......... + //...distance over voxel-hole line from voxel and hole to detection plane.......... - float dvd_l = (wmh.prj.rad - v->y1) / uy1; + float dvd_l = (wmh.prj.rad - v->y1) / uy1; - l->x1d_l = v->x1 + dvd_l * ux1; - l->z1d_l = v->z + dvd_l * uz1; + l->x1d_l = v->x1 + dvd_l * ux1; + l->z1d_l = v->z + dvd_l * uz1; - //...shadow of the hole .......... + //...shadow of the hole .......... - l->hsxcm_d = h->dxcm * dvd_l / dvh_l; - l->hszcm_d = h->dzcm * dvd_l / dvh_l; - l->hsxcm_d_d2 = l->hsxcm_d / (float)2.; - l->hszcm_d_d2 = l->hszcm_d / (float)2.; + l->hsxcm_d = h->dxcm * dvd_l / dvh_l; + l->hszcm_d = h->dzcm * dvd_l / dvh_l; + l->hsxcm_d_d2 = l->hsxcm_d / (float)2.; + l->hszcm_d_d2 = l->hszcm_d / (float)2.; - //... values at detection + crystal distance ................................ + //... values at detection + crystal distance ................................ - if (wmh.do_depth) - { + if (wmh.do_depth) + { - float dvdc_l = (wmh.prj.radc - v->y1) / uy1; + float dvdc_l = (wmh.prj.radc - v->y1) / uy1; - l->hsxcm_dc = h->dxcm * dvdc_l / dvh_l; - l->hszcm_dc = h->dzcm * dvdc_l / dvh_l; + l->hsxcm_dc = h->dxcm * dvdc_l / dvh_l; + l->hszcm_dc = h->dzcm * dvdc_l / dvh_l; - l->hsxcm_dc_d2 = l->hsxcm_d / (float)2.; - l->hszcm_dc_d2 = l->hszcm_d / (float)2.; + l->hsxcm_dc_d2 = l->hsxcm_d / (float)2.; + l->hszcm_dc_d2 = l->hszcm_d / (float)2.; - l->x1dc_l = v->x1 + dvdc_l * ux1; - l->z1dc_l = v->z + dvdc_l * uz1; + l->x1dc_l = v->x1 + dvdc_l * ux1; + l->z1dc_l = v->z + dvdc_l * uz1; + } + + //... effectiveness ...................................................... + + l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); } + else + { + // channel-edge margins + + float mxdx = h->x1 + h->dxcm / (float)2.; + float mxsx = h->x1 - h->dxcm / (float)2.; + float mzdx = h->z1 + h->dzcm / (float)2.; + float mzsx = h->z1 - h->dzcm / (float)2.; + float myup = h->y1 + h->dycm / (float)2.; + float mydn = h->y1 - h->dycm / (float)2.; + + // pinhole shadow for the region projecting inside the channel + if (v->x1 >= mxsx && v->x1 <= mxdx && v->z >= mzsx && v->z <= mzdx) + { + + //...vector voxel-hole, angles and distances............................... + + float ux1 = h->x1 - v->x1; + float uy1 = myup - v->y1; + float uz1 = h->z1 - v->z; + + if (uy1 <= EPSILON) error_weight3d(88, ""); + + //...vector voxel-hole and distances............................... + + float dxyz_2 = ux1 * ux1 + uy1 * uy1 + uz1 * uz1; + float dvh_l = sqrtf(dxyz_2); + + ux1 /= dvh_l; + uy1 /= dvh_l; + uz1 /= dvh_l; + + //...distance over voxel-hole line from voxel and hole to detection plane.......... + + float dvd_l = (wmh.prj.rad - v->y1) / uy1; + + l->x1d_l = v->x1 + dvd_l * ux1; + l->z1d_l = v->z + dvd_l * uz1; + + //...shadow of the hole .......... + + l->hsxcm_d = h->dxcm * dvd_l / dvh_l; + l->hszcm_d = h->dzcm * dvd_l / dvh_l; + l->hsxcm_d_d2 = l->hsxcm_d / (float)2.; + l->hszcm_d_d2 = l->hszcm_d / (float)2.; + //... values at detection + crystal distance ................................ + + if (wmh.do_depth) + { + + float dvdc_l = (wmh.prj.radc - v->y1) / uy1; + + l->hsxcm_dc = h->dxcm * dvdc_l / dvh_l; + l->hszcm_dc = h->dzcm * dvdc_l / dvh_l; + + l->hsxcm_dc_d2 = l->hsxcm_d / (float)2.; + l->hszcm_dc_d2 = l->hszcm_d / (float)2.; + + l->x1dc_l = v->x1 + dvdc_l * ux1; + l->z1dc_l = v->z + dvdc_l * uz1; + } + + //... effectiveness ...................................................... + + l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); + + } + else if ((v->x1 < mxsx || v->x1 > mxdx) && (v->z < mzsx || v->z > mzdx)) // x and z outside + { + float dxvch = (float)0.; + float dzvch = (float)0.; + + // voxel distance to the channel margin + if (v->x1 > 0) + { + dxvch = mxdx - v->x1; + } + else + { + dxvch = mxsx - v->x1; + } + if (v->z > 0) + { + dzvch = mzdx - v->z; + } + else + { + dzvch = mzsx - v->z; + } + + float dyvch = mydn - v->y1; + + // voxel to channel distance projected to the hole plane + + float dxh = dxvch * h->dycm / dyvch; + float dzh = dzvch * h->dycm / dyvch; - //... effectiveness ...................................................... + //...vector voxel-hole, angles and distances............................... - l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); + float ux1 = h->x1 + dxh - v->x1; + float uy1 = myup - v->y1; + float uz1 = h->z1 + dzh - v->z; + + if (uy1 <= EPSILON) + error_weight3d(88, ""); + + //...vector voxel-hole and distances............................... + + float dxyz_2 = ux1 * ux1 + uy1 * uy1 + uz1 * uz1; + float dvh_l = sqrtf(dxyz_2); + + ux1 /= dvh_l; + uy1 /= dvh_l; + uz1 /= dvh_l; + + //...distance over voxel-hole line from voxel and hole to detection plane.......... + + float dvd_l = (wmh.prj.rad - v->y1) / uy1; + + l->x1d_l = v->x1 + dvd_l * ux1; + l->z1d_l = v->z + dvd_l * uz1; + + //...shadow of the hole .......... + + l->hsxcm_d = (h->dxcm-abs(dxh)) * dvd_l / dvh_l; + l->hszcm_d = (h->dzcm-abs(dzh)) * dvd_l / dvh_l; + l->hsxcm_d_d2 = l->hsxcm_d / (float)2.; + l->hszcm_d_d2 = l->hszcm_d / (float)2.; + //... values at detection + crystal distance ................................ + + if (wmh.do_depth) + { + + float dvdc_l = (wmh.prj.radc - v->y1) / uy1; + + l->hsxcm_dc = (h->dxcm - abs(dxh)) * dvdc_l / dvh_l; + l->hszcm_dc = (h->dzcm - abs(dzh)) * dvdc_l / dvh_l; + + l->hsxcm_dc_d2 = l->hsxcm_d / (float)2.; + l->hszcm_dc_d2 = l->hszcm_d / (float)2.; + + l->x1dc_l = v->x1 + dvdc_l * ux1; + l->z1dc_l = v->z + dvdc_l * uz1; + } + + //... effectiveness ...................................................... + + l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); + } + else if ((v->x1 >= mxsx && v->x1 <= mxdx) && (v->z < mzsx || v->z > mzdx)) // x inside and z outside + { + float dzvch = (float)0.; + + if (v->z > 0) + { + dzvch = mzdx - v->z; + } + else + { + dzvch = mzsx - v->z; + } + + float dyvch = mydn - v->y1; + + // voxel to channel distance projected to the hole plane + + float dzh = dzvch * h->dycm / dyvch; + + //...vector voxel-hole, angles and distances............................... + + float ux1 = h->x1 - v->x1; + float uy1 = myup - v->y1; + float uz1 = h->z1 + dzh - v->z; + + if (uy1 <= EPSILON) + error_weight3d(88, ""); + + //...vector voxel-hole and distances............................... + + float dxyz_2 = ux1 * ux1 + uy1 * uy1 + uz1 * uz1; + float dvh_l = sqrtf(dxyz_2); + + ux1 /= dvh_l; + uy1 /= dvh_l; + uz1 /= dvh_l; + + //...distance over voxel-hole line from voxel and hole to detection plane.......... + + float dvd_l = (wmh.prj.rad - v->y1) / uy1; + + l->x1d_l = v->x1 + dvd_l * ux1; + l->z1d_l = v->z + dvd_l * uz1; + + //...shadow of the hole .......... + + l->hsxcm_d = h->dxcm * dvd_l / dvh_l; + l->hszcm_d = (h->dzcm - abs(dzh)) * dvd_l / dvh_l; + l->hsxcm_d_d2 = l->hsxcm_d / (float)2.; + l->hszcm_d_d2 = l->hszcm_d / (float)2.; + //... values at detection + crystal distance ................................ + + if (wmh.do_depth) + { + + float dvdc_l = (wmh.prj.radc - v->y1) / uy1; + + l->hsxcm_dc = h->dxcm * dvdc_l / dvh_l; + l->hszcm_dc = (h->dzcm - abs(dzh)) * dvdc_l / dvh_l; + + l->hsxcm_dc_d2 = l->hsxcm_d / (float)2.; + l->hszcm_dc_d2 = l->hszcm_d / (float)2.; + + l->x1dc_l = v->x1 + dvdc_l * ux1; + l->z1dc_l = v->z + dvdc_l * uz1; + } + + //... effectiveness ...................................................... + + l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); + } + else if ((v->x1 < mxsx || v->x1 > mxdx) && (v->z >= mzsx && v->z <= mzdx)) // x outside and z inside + { + float dxvch = (float)0.; + + if (v->x1 > 0) + { + dxvch = mxdx - v->x1; + } + else + { + dxvch = mxsx - v->x1; + } + + float dyvch = mydn - v->y1; + + // voxel to channel distance projected to the hole plane + + float dxh = dxvch * h->dycm / dyvch; + + //...vector voxel-hole, angles and distances............................... + + float ux1 = h->x1 + dxh - v->x1; + float uy1 = myup - v->y1; + float uz1 = h->z1 - v->z; + + if (uy1 <= EPSILON) + error_weight3d(88, ""); + + //...vector voxel-hole and distances............................... + + float dxyz_2 = ux1 * ux1 + uy1 * uy1 + uz1 * uz1; + float dvh_l = sqrtf(dxyz_2); + + ux1 /= dvh_l; + uy1 /= dvh_l; + uz1 /= dvh_l; + + //...distance over voxel-hole line from voxel and hole to detection plane.......... + + float dvd_l = (wmh.prj.rad - v->y1) / uy1; + + l->x1d_l = v->x1 + dvd_l * ux1; + l->z1d_l = v->z + dvd_l * uz1; + + //...shadow of the hole .......... + + l->hsxcm_d = (h->dxcm - abs(dxh)) * dvd_l / dvh_l; + l->hszcm_d = h->dzcm * dvd_l / dvh_l; + l->hsxcm_d_d2 = l->hsxcm_d / (float)2.; + l->hszcm_d_d2 = l->hszcm_d / (float)2.; + //... values at detection + crystal distance ................................ + + if (wmh.do_depth) + { + + float dvdc_l = (wmh.prj.radc - v->y1) / uy1; + + l->hsxcm_dc = (h->dxcm - abs(dxh)) * dvdc_l / dvh_l; + l->hszcm_dc = h->dzcm * dvdc_l / dvh_l; + + l->hsxcm_dc_d2 = l->hsxcm_d / (float)2.; + l->hszcm_dc_d2 = l->hszcm_d / (float)2.; + + l->x1dc_l = v->x1 + dvdc_l * ux1; + l->z1dc_l = v->z + dvdc_l * uz1; + } + + //... effectiveness ...................................................... + + l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); + } + + + } } //========================================================================== From a1b56cea12833b8098e8e0b792e3e53a992aeaaf Mon Sep 17 00:00:00 2001 From: TommasoFerri97 <79870251+TommasoFerri97@users.noreply.github.com> Date: Wed, 2 Sep 2026 13:33:58 +0200 Subject: [PATCH 2/3] Fix channel-edge projected centroid calculation --- src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx | 16 ++++++++-------- 1 file changed, 8 insertions(+), 8 deletions(-) diff --git a/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx b/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx index b05655c023..59edf4b7e4 100644 --- a/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx +++ b/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx @@ -602,7 +602,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const float dzvch = (float)0.; // voxel distance to the channel margin - if (v->x1 > 0) + if (v->x1 > h->x1) { dxvch = mxdx - v->x1; } @@ -610,7 +610,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const { dxvch = mxsx - v->x1; } - if (v->z > 0) + if (v->z > h->z1) { dzvch = mzdx - v->z; } @@ -628,9 +628,9 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const //...vector voxel-hole, angles and distances............................... - float ux1 = h->x1 + dxh - v->x1; + float ux1 = h->x1 + dxh / 2.F - v->x1; float uy1 = myup - v->y1; - float uz1 = h->z1 + dzh - v->z; + float uz1 = h->z1 + dzh / 2.F - v->z; if (uy1 <= EPSILON) error_weight3d(88, ""); @@ -682,7 +682,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const { float dzvch = (float)0.; - if (v->z > 0) + if (v->z > h->z1) { dzvch = mzdx - v->z; } @@ -701,7 +701,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const float ux1 = h->x1 - v->x1; float uy1 = myup - v->y1; - float uz1 = h->z1 + dzh - v->z; + float uz1 = h->z1 + dzh / 2.F - v->z; if (uy1 <= EPSILON) error_weight3d(88, ""); @@ -753,7 +753,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const { float dxvch = (float)0.; - if (v->x1 > 0) + if (v->x1 > h->x1) { dxvch = mxdx - v->x1; } @@ -770,7 +770,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const //...vector voxel-hole, angles and distances............................... - float ux1 = h->x1 + dxh - v->x1; + float ux1 = h->x1 + dxh / 2.F - v->x1; float uy1 = myup - v->y1; float uz1 = h->z1 - v->z; From b5dbbb097e00c3942464f1b5bb1061eea1f62633 Mon Sep 17 00:00:00 2001 From: TommasoFerri97 <79870251+TommasoFerri97@users.noreply.github.com> Date: Wed, 2 Sep 2026 16:00:44 +0200 Subject: [PATCH 3/3] Fixed channel-edge area fraction calculation --- .../PinholeSPECTUB_Weight3d.cxx | 44 +++++++++++++++++++ 1 file changed, 44 insertions(+) diff --git a/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx b/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx index 59edf4b7e4..a54a5ab567 100644 --- a/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx +++ b/src/recon_buildblock/PinholeSPECTUB_Weight3d.cxx @@ -49,6 +49,8 @@ static bool check_zang_par(const voxel_type* vox, const hole_type* h); // bool check_zang_obl( lor_type * l, voxel_type * vox, hole_type * h); +static float channel_area_fraction(const hole_type* h, const float dxh, const float dzh); + static void voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const wmh_mph_type& wmh); static void downsample_psf(const psf2d_type* psf_in, psf2d_type* psf_out, int factor, bool do_calc); @@ -467,6 +469,45 @@ check_zang_par(const voxel_type* v, const hole_type* h) return (ans); } +//========================================================================== +//=== channel_area_fraction ================================================ +//========================================================================== + +static float +channel_area_fraction(const hole_type* h, const float dxh, const float dzh) +{ + if (h->do_round) + { + // Relative displacement of two identical elliptical/circular apertures. + const float qx = dxh / h->dxcm; + const float qz = dzh / h->dzcm; + const float q = (std::sqrt)(qx * qx + qz * qz); + + if (q >= 1.F) + return 0.F; + + constexpr float pi = 3.14159265358979323846F; + + const float radicand = 1.F - q * q; + + return (2.F / pi) + * ((std::acos)(q) + - q * (std::sqrt)(radicand > 0.F ? radicand : 0.F)); + } + + // Rectangular aperture. + const float fx_raw + = 1.F - (std::abs)(dxh) / h->dxcm; + + const float fz_raw + = 1.F - (std::abs)(dzh) / h->dzcm; + + const float fx = fx_raw > 0.F ? fx_raw : 0.F; + const float fz = fz_raw > 0.F ? fz_raw : 0.F; + + return fx * fz; +} + //========================================================================== //=== voxel_projection ===================================================== //========================================================================== @@ -677,6 +718,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const //... effectiveness ...................................................... l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); + l->eff *= channel_area_fraction(h, dxh, dzh); } else if ((v->x1 >= mxsx && v->x1 <= mxdx) && (v->z < mzsx || v->z > mzdx)) // x inside and z outside { @@ -748,6 +790,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const //... effectiveness ...................................................... l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); + l->eff *= channel_area_fraction(h, 0.F, dzh); } else if ((v->x1 < mxsx || v->x1 > mxdx) && (v->z >= mzsx && v->z <= mzdx)) // x outside and z inside { @@ -819,6 +862,7 @@ voxel_projection_mph(lor_type* l, const voxel_type* v, const hole_type* h, const //... effectiveness ...................................................... l->eff = wmh.mndvh2 / dxyz_2 * fabsf(uy1); + l->eff *= channel_area_fraction(h, dxh, 0.F); }