Problème programme Image médicale PET/CT
Bonjour,
J'aurais besoin d'un peut d'aide car je suis devant un problème.
Je suis entrain de modifier une application pour l'imagerie médicale.
elle permet en temps normal de calculer des facteurs de correction de l'attenuation a partir d'image obtenu a partir d'un scanner X.
La partie du programme qui me pose un problème est la projection qui permet d'obtenir a partir d'une image CT (scanner X), le sinogramme correspondant, ce code (écrit par quelqu'un d'autre) ne fonctionne pas et j'éssaie depuis plus d'une semaine de la debugger.
Code:
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
|
void CDAttnCT::proj_lookup(void)
{
float l=float(max_t-1)*0.5f;
for (int iphi=0; iphi<max_phi*6; iphi+=4)
{
float phi=float(iphi/4)*M_PI/float(max_phi*1.5);
geo_lookup[iphi]=sin(phi);
geo_lookup[(iphi)+1]=cos(phi);
geo_lookup[(iphi)+2]=l*(1.f+geo_lookup[iphi]-geo_lookup[iphi+1]);
geo_lookup[(iphi)+3]=-l*(geo_lookup[iphi]+geo_lookup[iphi+1]);
}
float halfmaxt2=l*l;
for (int i=1; i<max_t-1; i++) {
geo_factor[i]=1.f/sqrt(halfmaxt2-(float(i)-l)*(float(i)-l));
}
}
projection directe vers les sinogrammes
void CDAttnCT::proj_dir(void)
{
float margin=float(max_t-1);//max_t =255
float radius=ring_radius/bin_size;//ring_radius =87.3 et bin_size c'est 0.4
float radius2=radius*radius;
float halfmax_t=margin*0.5f;
float imagevalue;
int slice;
for (slice=0; slice<no_img_slice; slice+=1) //only direct image planes
{//no_img_slice c'est 95
int ind= (slice)*max_t*max_phi; // for direct sinos
int indimage=slice*max_t2;
//max_phi=120
for (int phi=0; phi<max_phi*1.5; phi++)
{
float sinphi=geo_lookup[(phi)<<2];
float cosphi=geo_lookup[((phi)<<2)+1];
float f1=geo_lookup[((phi)<<2)+2];
int indx=ind+phi*max_t;
//for (int yi=low_l; yi<high_l; yi++)
for (int yi=0; yi<max_t; yi++)
{//max_t=255
float yyt=f1-float(yi)*sinphi;
int indimage2=yi*max_t+indimage;
//for (int xi=low_l; xi<high_l; xi++)
for (int xi=0; xi<max_t; xi++)
// ignore negative part
if ((imagevalue=image[indimage2+xi])!=INVALID_V2 && imagevalue>0.f)
{
float t=yyt+float(xi)*cosphi;
if (t<0.f || t>margin)
continue;
int tbi=int(t);
float tl=t-float(tbi);
float th=1.f-tl;
dir_reproj[indx+tbi]+=imagevalue*th;
dir_reproj[indx+tbi+1]+=imagevalue*tl;
}
}
}
}
} |
A la sortie, le sinogramme est completement faux.
Est ce que quelqu'un pourrait m'aider et me dire si quelque chose ne tient pas la route dans ce code?
La premiere partie calcule apparemment des coefficient géométriques
et la deuxiéme calcule les projections.
Merci d'avance
