Bonsoir,

Je suis étudiante, et je travail actuellement sur un projet en C++ : la reconstruction d'une image fantôme.
Je dois effectuer une reconstruction tomographique en effectuant une rétroprojection des coefficients d'atténuations moyens mesurés
Le code à modifier est le suivant : (programme principal main qui fait appel à des classes ROOT)

Code C++ : Sélectionner tout - Visualiser dans une fenêtre à part
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
#include <iomanip>
#include <time.h>
#include <signal.h>
#include <math.h>
 
#include <TMath.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>
#include <TF2.h>
#include <TH1F.h>
#include <TH1D.h>
#include <TH2D.h>
#include <TH3D.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <THStack.h>
#include <TLegend.h>
#include <TROOT.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TTree.h>
#include <TArrayD.h>
#include <TMatrixD.h>
#include <TVectorD.h>
#include <TBranch.h>
#include <TColor.h>
#include <TRandom3.h>
#include <TFile.h>
#include <TApplication.h>
#include <TProfile2D.h>
#include <TRandom3.h>
 
using namespace std;
 
 
int main(int argc, char **argv)
{
 
  //#########################################################
  // root handling
  // ./main -b : batch mode (no display of the figures)
  //#########################################################
  TApplication theApp("App", &argc, argv); 
 
  // Factor of statistics reduction (to speed up the program execution)
  Double_t StatFactor = 1;     
 
  //root colz style (smooth 2D palette)
  const Int_t NRGBs = 5;
  const Int_t NCont = 255;
  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
  gStyle->SetNumberContours(NCont);
 
  //General settings
  static const Int_t radius = 100; //phantom (cylinder) radius in mm
  Double_t I0 = 500.; //photons number per bin for each projection 
  Double_t pixelSize = 1.; //in mm
  static const Int_t Nproj = 180; //projection number
  Double_t theta[Nproj] = {0.};
  Double_t AngularStep = 1.; // Angular step bewteen the radiography
 
  for (int div=0; div<Nproj; div++) theta[div] = AngularStep*div*TMath::Pi()/180.; //RunID = projection angle
 
//   Int_t iter = 3; //iteration number
  static const Int_t idiam = 2*radius;
 
// Histograms
 
//   TH2F *hEnergyZ = new TH2F("hEnergyZ","hEnergyZ",200,-100,100,100,0,100);
  TH2F *hYZ = new TH2F("hYZ","hYZ",200,-100,100,21,-10,10);	
 
 
  //CT Tables x=angle, y=position on projection axis
  Double_t TabCT[Nproj][idiam] = {{0.,0.}};
 
  //Get tree results (Gate output)
 
  TFile *f = new TFile("RootFile/PhaseSpace.root","read");
  TTree *tree      = (TTree *) f->Get("PhaseSpace");
 
  //declare tree variables
  Float_t Y, Z, Ekine; // position in mm and Ekine in eV
  Int_t RunID;
 
  //Get tree variables
  tree->SetBranchAddress("Z",&Z);
  tree->SetBranchAddress("Y",&Y);
  tree->SetBranchAddress("Ekine",&Ekine);
  tree->SetBranchAddress("RunID", &RunID );
 
  //loop over tree entries
  Int_t nentries = (Int_t)tree->GetEntries();
  printf("Number of entries = %d\n",nentries);
 
  cout << "Number of read entries: " << nentries/StatFactor << endl;
 
  for (Int_t n=0; n < nentries/StatFactor; n++) 
  {							
    //Get tree entry n
    tree->GetEntry(n);
 
    // Energy E(keV) and position (mm) on the detector
//     hEnergyZ->Fill(Z,Ekine*1000.,1.);
    hYZ->Fill(Y,Z,1.);		
  }
 
  gStyle->SetOptStat(0000);
 
  TCanvas *c0 = new TCanvas("hYZ","hYZ",500,500);
  hYZ->Draw("colz");
  hYZ->GetXaxis()->SetTitle("Position Y (mm)");
  hYZ->GetYaxis()->SetTitle("Position Z (mm)");
  c0->SetLogz();
 
//   TCanvas *c1 = new TCanvas("hEnergyZ","hEnergyZ",500,500);
//   hEnergyZ->Draw("colz");
//   hEnergyZ->GetXaxis()->SetTitle("Position Z (mm)");
//   hEnergyZ->GetYaxis()->SetTitle("Energie (keV)");
//   c1->SetLogz();
 
  // Save the figures
  c0-> SaveAs("Figures/hZY.pdf");
//   c1-> SaveAs("Figures/hEnergyZ.pdf");
 
  // Display of the figures if "interactive" mode (not batch)
  if(!gROOT->IsBatch()) theApp.Run();
 
  return 0;
 
}

Je dois déterminer la rétroprojection filtrée de la première radiographique, puis généraliser le code pour l'ensemble des radiographies: je ne sais pas comment m'y prendre.

Merci pour votre aide