first commit

.gitignore

+build
\ No newline at end of file

.vscode/c_cpp_properties.json

+{
+    "configurations": [
+        {
+            "name": "Linux",
+            "includePath": [
+                "${workspaceFolder}/**",
+                "/home/ming/root/include",
+                "/home/ming/root/lib"
+            ],
+            "defines": [],
+            "compilerPath": "/usr/bin/gcc",
+            "cStandard": "gnu17",
+            "cppStandard": "gnu++14",
+            "intelliSenseMode": "linux-gcc-x64",
+            "configurationProvider": "ms-vscode.cmake-tools"
+        }
+    ],
+    "version": 4
+}
\ No newline at end of file

CMakeLists.txt

+cmake_minimum_required(VERSION 3.16.3)
+
+project(Imager CXX)
+
+set(CMAKE_CXX_STANDARD 11)
+set(CMAKE_POSITION_INDEPENDENT_CODE ON)
+
+# find_package(ROOT 6.22 CONFIG REQUIRED)
+find_package(ROOT 6.22 CONFIG REQUIRED COMPONENTS Gui)
+
+include_directories(${CMAKE_CURRENT_SOURCE_DIR}/Headers)
+
+enable_testing()
+# add_subdirectory(googletest)
+add_subdirectory(Sources)
+# add_subdirectory(unitest)
+# add_subdirectory(bin)
\ No newline at end of file

Headers/draw.h

+#pragma once
+
+#include <iostream>
+
+#include "TCanvas.h"
+#include "TH2F.h"
+#include "TGraph.h"
+#include "TStyle.h"
+#include "TEllipse.h"
+
+#include "../Headers/setup.h"
+
+int Backprojection(const Setup& config, const std::vector<Cone>& cones, 
+                   std::vector<std::vector<double>>& image);
+
+int drawImage(const Setup& config, const std::vector<std::vector<double>>& image);
\ No newline at end of file

Headers/setup.h

+#pragma once
+
+#include <vector>
+#include <math.h>
+#include <stdexcept>
+
+class Vector3D
+{
+public:
+    double X=0;
+    double Y=0;
+    double Z=0;
+    Vector3D():Vector3D(0, 0, 0) {}
+    Vector3D(const double x, const double y, const double z): X(x), Y(y), Z(z) {}
+
+    Vector3D operator+(const Vector3D& r) const {
+        return Vector3D(this->X + r.X, this->Y + r.Y, this->Z + r.Z);
+    }
+    Vector3D operator-(const Vector3D& r) const {
+        return Vector3D(this->X - r.X, this->Y - r.Y, this->Z - r.Z);
+    }
+
+    double operator*(const Vector3D& r) const {
+        return this->X * r.X + this->Y * r.Y + this->Z * r.Z;
+    }
+
+};
+
+
+inline double getCosAngle(const Vector3D& l, const Vector3D& r) {
+    return l*r / std::sqrt((l*l) * (r*r));
+}
+
+class Setup
+{
+private:
+    /* data */
+public:
+    double pitchX = 11;
+    int      numX =  7;
+    double pitchY = 13;
+    int      numY =  4;
+    double pitchZ = 35;
+    int      numZ =  2;
+    // pixel coordinates, mm
+    std::vector<Vector3D> pixels;
+    int getPixelPos();
+    // radius of projection sphere,  mm
+    double R = 500;
+    Setup(): Setup(500) {}
+    Setup(const double r):R(r) {
+        getPixelPos(); 
+        pixelateSphere();
+        randSource();
+    }
+
+    // image setup
+    // Altitude angle, rad
+    const int thetaBins = 179;
+    const int phiBins = 360;
+    const double dtheta = 179 / double(thetaBins) * M_PI / 180;
+    const double thetaMin = -89 * M_PI / 180;
+    const double thetaMax = 89 * M_PI / 180;
+    // Azimuthal angle, rad
+    const double dphi = 2 * M_PI /double(phiBins);
+    const double phiMin = -M_PI + dphi / 2;
+    const double phiMax = M_PI - dphi / 2;
+    // pixelate the projection sphere
+    std::vector<std::vector<Vector3D>> xbs;
+    int pixelateSphere();
+
+    // simulation setup
+    const int nCones= 1000;
+    double trueTheta;
+    double truePhi;
+    Vector3D trueSource;
+    int randSource();
+};
+
+class Cone
+{
+public:
+    // interaction site 1 is the apex
+    Vector3D apex;
+    // site 2 - site 1 is the axis
+    Vector3D axis;
+    // half angle = abs (scattering angle)
+    // for neutron, 0 < half angle < pi /2
+    // for gamma,   0 < half angle < pi
+    double cosHalfAngle;
+    Cone(const Vector3D point, const Vector3D line, const double t): apex(point), axis(line), cosHalfAngle(t) {}
+};
+
+

Headers/simulation.h

+#pragma once
+
+#include <cstdlib>
+#include <ctime>
+
+#include  "setup.h"
+
+Cone createCone(const Setup& config);
\ No newline at end of file

Sources/CMakeLists.txt

+add_library(setup setup.cpp)
+
+add_library(simulation simulation.cpp)
+target_link_libraries(simulation PUBLIC setup)
+
+add_library(draw draw.cpp)
+target_link_libraries(draw PUBLIC setup ROOT::Gpad)
+
+add_executable(main main.cpp)
+target_link_libraries(main PUBLIC setup simulation draw ROOT::Gpad ROOT::Gui)
+# set_target_properties(main PROPERTIES RUNTIME_OUTPUT_DIRECTORY "${CMAKE_CURRENT_SOURCE_DIR}")
\ No newline at end of file

Sources/draw.cpp

+#include "../Headers/draw.h"
+
+int Backprojection(const Setup& config, const std::vector<Cone>& cones, 
+                   std::vector<std::vector<double>>& image){
+    // project cones onto the spherical surface
+    std::cout << "Projecting cones onto the designated spherical surface..."<<std::endl;
+    double progress=0.0;
+    u_int currentConeNum(0);
+    u_int totalConeNum(cones.size());
+    double alpha(0);
+    double beta(0);
+    for (const Cone& event : cones)
+    {
+        // update progress bar
+        int barWidth = 70;
+        currentConeNum++;
+        if (currentConeNum % 10 == 0)
+        {
+            progress = double(currentConeNum) / double(totalConeNum);
+            std::cout << "[";
+            int pos = barWidth * progress;
+            for (int i = 0; i < barWidth; ++i) {
+                if (i < pos) std::cout << "=";
+                else if (i == pos) std::cout << ">";
+                else std::cout << " ";
+            }
+            std::cout << "] " << int(progress * 100.0) << " %\r";
+            std::cout.flush();
+        }
+
+        // update the image
+        //TODO: add the uncertainty
+        alpha = event.cosHalfAngle;
+        for (int i = 0; i < config.thetaBins; i++)
+        {
+            for (int j = 0; j < config.phiBins; j++)
+            {
+                beta = getCosAngle(event.apex - config.xbs[i][j], event.axis);
+                image[i][j] += std::exp(-(beta-alpha)*(beta-alpha) / 0.002);
+            }
+        }
+    }
+    return 0;
+}
+
+int drawImage(const Setup& config, const std::vector<std::vector<double>>& image){
+    TH2D *histo = new TH2D("ROI", " ; Azimuth; Elevation", 
+                          config.phiBins, -180, 180, 
+                          config.thetaBins, -89.5, 89.5);
+    TCanvas *canvas = new TCanvas("THCanvas","Sinocanvas", 1000, 500);
+    for (int i = 1; i <= config.phiBins; i++)
+    {
+        for (int j = 1; j < config.thetaBins; j++)
+        {
+            histo->SetBinContent(i, j, image[j][i]);
+        }
+    }
+    gStyle->SetOptStat(0);
+    histo->GetZaxis()->SetLabelSize(0.02);
+    histo->Draw("z aitoff");
+
+    // grid
+    float conv=M_PI/180; // I am also aware of TMath::DegToRad() and TMath::Pi() which could be used there...
+    float la, lo, x, y, z;
+    
+    const int Nl = 19; // Number of drawn latitudes
+    const int NL = 19; // Number of drawn longitudes
+    int M  = 30;
+
+    TGraph  *latitudes[Nl];
+    TGraph  *longitudes[NL];
+
+    for (int j=0;j<Nl;++j) {
+        latitudes[j]=new TGraph();
+        la = -90+180/(Nl-1)*j;
+        for (int i=0;i<M+1;++i) {
+            lo = -180+360/M*i;
+            z  = sqrt(1+cos(la*conv)*cos(lo*conv/2));
+            x  = 180*cos(la*conv)*sin(lo*conv/2)/z;
+            y  = 90*sin(la*conv)/z;
+            latitudes[j]->SetPoint(i,x,y);
+        }
+    } 
+    
+    for (int j=0;j<NL;++j) {
+        longitudes[j]=new TGraph();
+        lo = -180+360/(NL-1)*j;
+        for (int i=0;i<M+1;++i) {
+            la = -90+180/M*i;
+            z  = sqrt(1+cos(la*conv)*cos(lo*conv/2));
+            x  = 180*cos(la*conv)*sin(lo*conv/2)/z;
+            y  = 90*sin(la*conv)/z;
+            longitudes[j]->SetPoint(i,x,y);
+        }
+    }
+    // Draw the grid	
+    TPad *pad2 = new TPad("pad2","",0,0,1,1);
+    pad2->SetFillStyle(4000);
+    pad2->SetFillColor(0);
+    pad2->SetBorderSize(0);
+    pad2->SetFrameBorderMode(0);
+    pad2->SetFrameLineColor(0); 
+    pad2->SetFrameBorderMode(0);
+    pad2->Draw();
+    pad2->cd();
+    Double_t ymin = -89.5;
+    Double_t ymax = 89.5;
+    Double_t dy = (ymax-ymin)/0.8; //10 per cent margins top and bottom
+    Double_t xmin = -180;
+    Double_t xmax = 180;
+    Double_t dx = (xmax-xmin)/0.8; //10 per cent margins left and right
+
+    pad2->Range(xmin-0.1*dx,ymin-0.1*dy,xmax+0.1*dx,ymax+0.1*dy);
+
+    for (int j=0;j<Nl;++j) latitudes[j]->Draw("l");
+    for (int j=0;j<NL;++j) longitudes[j]->Draw("l");
+    // // // canvas1->cd(4)->SetLogz();
+
+    // // // draw source 1
+    // // TEllipse* el1 = new TEllipse(86.5276, 0, 5., 5.);
+    // // el1->SetFillColor(0);
+    // // el1->SetFillStyle(0);
+    // // el1->SetLineColor(4);
+    // // el1->Draw("SAME");
+
+    // draw source 2
+    double sigma = 5;
+    TEllipse* el2 = new TEllipse(config.truePhi / conv, config.trueTheta / conv, sigma, sigma);
+    el2->SetFillColor(0);
+    el2->SetFillStyle(0);
+    el2->SetLineColor(4);
+    el2->Draw("SAME");
+    canvas->Draw();
+
+    return 0;
+}
\ No newline at end of file

Sources/main.cpp

+#include <vector>
+#include <iostream>
+#include <cstdlib>
+#include <ctime>
+
+#include "TApplication.h"
+
+#include "../Headers/setup.h"
+#include "../Headers/simulation.h"
+#include "../Headers/draw.h"
+
+int main(int argc, char** argv) {
+    // init
+    const Setup config(10);
+    // simulate some cones
+    std::vector<Cone> cones;
+    std::srand(std::time(nullptr));
+    // for (int i = 0; i < 10; i++)
+    for (int i = 0; i < config.nCones; i++)
+    {
+        cones.push_back(createCone(config));
+    }
+    // generate a projection image
+    // init image
+    std::vector<std::vector<double>> image(config.thetaBins, std::vector<double>(config.phiBins, 0));
+    // draw    
+    TApplication *myapp = new TApplication("myApp", 0, 0);
+    Backprojection(config, cones, image);
+    drawImage(config, image);
+    myapp->Run();
+    return 0;
+}
\ No newline at end of file

Sources/setup.cpp

+#include "../Headers/setup.h"
+
+int Setup::getPixelPos()
+{
+    // get sipm pixel pos
+    for (int i = 0; i < numZ; i++)
+    {
+        double zi = i * pitchZ - (numZ - 1) * pitchZ / 2.0;
+        for (int j = 0; j < numX; j++)
+        {
+            double xj = j * pitchX - (numX - 1) * pitchX / 2.0;
+            for (int k = 0; k < numY; k++)
+            {
+                double yk = k * pitchY - (numY - 1) * pitchY / 2.0;
+                pixels.push_back(Vector3D(xj, yk, zi));
+            }
+        }
+    }
+    return 0;
+}
+
+int Setup::pixelateSphere() {
+    xbs.resize(thetaBins, std::vector<Vector3D>(phiBins, Vector3D(0, 0, 0)));
+    double phi = phiMin;
+    double theta = thetaMin;
+    for (int i = 0; i < thetaBins; i++)
+    {
+        for (int j = 0; j < phiBins; j++)
+        {
+            xbs[i][j].X = R * std::cos(theta) * std::cos(phi);
+            xbs[i][j].Y = R * std::cos(theta) * std::sin(phi);
+            xbs[i][j].Z = R * std::sin(theta);
+            phi += dphi;
+        }
+        theta += dtheta;
+    }
+    
+    return 0;
+}
+
+int Setup::randSource(){
+    // random  theta and phi
+    trueTheta = float(std::rand()) / RAND_MAX * (thetaMax - thetaMin) + thetaMin;
+    truePhi = float(std::rand()) / RAND_MAX * (phiMax - phiMin) + phiMin;
+    trueSource = Vector3D(R*std::cos(trueTheta)*std::cos(truePhi), 
+                          R*std::cos(trueTheta)*std::sin(truePhi),
+                          R*std::sin(trueTheta));
+    return 0;
+}
\ No newline at end of file

Sources/simulation.cpp

+#include "../Headers/simulation.h"
+
+Cone createCone(const Setup& config){
+    // static Vector3D fakeSource(config.R, 0, 0);
+
+    // randomly pick up two pixels
+    int i = std::rand() % (config.pixels.size() / 2);
+    int j = -1;
+    while (true)
+    {
+        j = std::rand() % (config.pixels.size() / 2);
+        if (j != i)
+        {
+            break;
+        }
+    }
+    Vector3D site1 = config.pixels[i];
+    Vector3D site2 = config.pixels[j];
+    // randomly select z coordinates of two interations
+    site1.Z = (2 * float(std::rand()) / RAND_MAX - 1) * config.pitchZ / 2.0;
+    site2.Z = (2 * float(std::rand()) / RAND_MAX - 1) * config.pitchZ / 2.0;
+    // calculate angle
+    Vector3D axis = site2 - site1;
+    double cosAngle = getCosAngle(site1 - config.trueSource, axis);
+    return Cone(site1, axis, cosAngle);
+}
\ No newline at end of file