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 : 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
 
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