GPGPU labor XII Tomogrfis rekonstrukci Kezdeti teendk Tantrgy

GPGPU labor XII. Tomográfiás rekonstrukció

Kezdeti teendők • Tantárgy honlapja, Monte Carlo szimuláció • A labor kiindulási alapjának letöltése (lab 12_base. zip), kitömörítés a D: GPGPU könyvtárba • D: GPGPUlabslab 12_opencllab 12_openc l. sln indítása • Project tulajdonságai – Configuration Properties – Debugging – Working Directory = $(Project. Dir). . bin

„Keretrendszer” • A „referencia mérés” legyártását és az EM iteráció vázát tartalmazza, üres projekciós kernelekkel • A „referencia mérés” modul normálissal párhuzamos LOR-jait (ld. később) képként kimenti: out_000. tga • „i”-re lépteti az iterációt – Képként elmenti a modul normálissal párhuzamos LORokat (ld. később): out_{iterációszám}. tga – Frissíti a volume-ot (a kontrasztanyag eloszlásának aktuális becslését) • Volume megjelenítő kóddal kiegészítve, az előző gyakorlat alapján – Hogy lássuk, mi az aktuális volume becslés

Modul geometria HOST: struct module { cl_float 4 }; origin; axial; trans. Axial; n; Device: struct module { float 4 origin; float 4 axial; float 4 trans. Axial; float 4 n; }; a 1 a 2 1 o 1 n 2 t 1 n 1 0 t 0 a 0 o 0 n 0 y o 2 z t 2 2 (axiális) x n 3 3 a 3 t 3 3

Modul geometria const int n. Modules = 4; const int coincidence = 1; const int n. Pairs = n. Modules*coincidence / 2; const int n. Axial = 32; const int n. Trans. Axial = 32; const float detector. Area = 2. 0 f*2. 0 f / n. Axial / n. Trans. Axial; module. Buffer. CPU[n. Modules]; cl_mem module. Buffer. GPU; void rotate 90 xy( cl_float 4 &vector ) { std: : swap( vector. s[0], vector. s[1] ); vector. s[0] *= -1. 0 f; } void init. Simulation(){ // … cl_float 4 origin, axial, trans. Axial, n; origin. s[3] = axial. s[3] = trans. Axial. s[3] = n. s[3] = 0. 0 f; // (x, x, x, 0) origin. s[0] = origin. s[1] = origin. s[2] = -1. 0 f; // (-1, -1) axial. s[2] = 2. 0 f; axial. s[0] = axial. s[1] = 0. 0 f; // (0, 0, 2) trans. Axial. s[0] = 2. 0 f; trans. Axial. s[1] = trans. Axial. s[2] = 0. 0 f; // (2, 0, 0) n. s[1] = 1. 0 f; n. s[0] = n. s[2] = 0. 0 f; // (0, 1, 0) detector. Area = 2. 0 f*2. 0 f / n. Axial / n. Trans. Axial; for ( int i = 0; i < n. Modules; ++i ) { module. Buffer. CPU[i]. origin = origin; module. Buffer. CPU[i]. axial = axial; module. Buffer. CPU[i]. trans. Axial = trans. Axial; module. Buffer. CPU[i]. n = n; rotate 90 xy(origin); rotate 90 xy(axial); rotate 90 xy(trans. Axial); rotate 90 xy(n); } CL_SAFE_CALL( cl. Enqueue. Write. Buffer(commands, module. Buffer. GPU, CL_TRUE, 0, sizeof(module)*n. Modules, module. Buffer. CPU, 0, NULL)); }

Ellenőrzés std: : ostream& operator<< ( std: : ostream& os , const cl_float 4& v ) { os << "(" << v. s[0] << ", " << v. s[1] << ", " << v. s[2] << ", " << v. s[3] << ")"; return os; } // pl. init-be: for ( unsigned int i { std: : cout << i << << ", t=" << ", o=" << ", n=" } 0: 1: 2: 3: = 0; i < n. Modules; ++i ) ": << << << a=" << module. Buffer. CPU[i]. axial module. Buffer. CPU[i]. trans. Axial module. Buffer. CPU[i]. origin module. Buffer. CPU[i]. n << std: : endl ; a=(0, 0, 2, 0), t=(2, 0, 0, 0), o=(-1, -1, 0), n=(0, 1, 0, 0) a=(-0, 0, 2, 0), t=(-0, 2, 0, 0), o=(1, -1, 0), n=(-1, 0, 0, 0) a=(-0, 2, 0), t=(-2, -0, 0, 0), o=(1, 1, -1, 0), n=(-0, -1, 0, 0) a=(0, -0, 2, 0), t=(0, -2, 0, 0), o=(-1, 1, -1, 0), n=(1, -0, 0, 0) a 1 a 2 1 o 1 n 2 t 1 n 1 0 t 0 a 0 o 0 n 0 y o 2 z t 2 2 (axiális) x n 3 3 a 3 t 3 3

Koincidencia • A host oldalon nincs rá szükség • A kernelekben is csak a párindex → panelindex irányra van szükség: int 2 get. Modules(int pair){ switch ( pair ) { case 0: return (int 2)(0, 2); case 1: return (int 2)(1, 3); }; return (int 2)(0, 0); }

Volume fizikai elhelyezkedése cl_float 4 volume. Min; cl_float 4 volume. Max; void init. Simulation(){ // … volume. Min. s[0] = volume. Min. s[1] = volume. Min. s[2] = -0. 5 f; volume. Max. s[0] = volume. Max. s[1] = volume. Max. s[2] = 0. 5 f; volume. Min. s[3] = volume. Max. s[3] = 0. 0 f; //… } volume. Max x volume. Min y z

Adatok // LOR const int lor. Number. Pair = n. Axial*n. Trans. Axial; const int lor. Number = lor. Number. Pair*n. Pairs; cl_mem estimated. Lor. Buffer. GPU; cl_mem measured. Lor. Buffer. GPU; float estimated. Lor. Buffer. CPU[lor. Number]; float measured. Lor. Buffer. CPU[lor. Number]; void init. Simulation(){ // … estimated. Lor. Buffer. GPU = cl. Create. Buffer(context, CL_MEM_READ_WRITE, sizeof(float) * lor. Number, NULL); measured. Lor. Buffer. GPU = cl. Create. Buffer(context, CL_MEM_READ_WRITE, sizeof(float) * lor. Number, NULL); // … } // Volume const int resolution = 32; cl_mem volume. Buffer. GPU; float* volume. Data; void init. Simulation(){ // … volume. Data = new float[resolution * resolution]; for ( unsigned int x = 0; x < resolution; ++x ) for ( unsigned int y = 0; y < resolution; ++y ) for ( unsigned int z = 0; z < resolution; ++z ) volume. Data[ x +resolution*y +resolution*z ] = 1. 0 f; volume. Buffer. GPU = cl. Create. Buffer(context, CL_MEM_READ_WRITE, sizeof(float) * resolution, NULL); CL_SAFE_CALL( cl. Enqueue. Write. Buffer(commands, volume. Buffer. GPU, CL_TRUE, 0, sizeof(float) * resolution, volume. Data, 0, NULL)); // … }

Vetítő operátorok • Előrevetítés • Visszavetítés

Előrevetítés: problématér felosztása // Host: void simulation. Step(){ CL_SAFE_CALL( cl. Set. Kernel. Arg(forward. Projection. Kernel, CL_SAFE_CALL( cl. Set. Kernel. Arg(forward. Projection. Kernel, 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, sizeof(int), &n. Axial) ); sizeof(int), &n. Trans. Axial) ); sizeof(int), &resolution) ); sizeof(float), &detector. Area) ); sizeof(cl_float 4), &volume. Min) ); sizeof(cl_float 4), &volume. Max) ); sizeof(cl_mem), &module. Buffer. GPU) ); sizeof(cl_mem), &measured. Lor. Buffer. GPU) ); sizeof(cl_mem), &estimated. Lor. Buffer. GPU) ); sizeof(cl_mem), &volume. Buffer. GPU) ); size_t forward. Projection. Size[3]; forward. Projection. Size[0] = n. Pairs; forward. Projection. Size[1] = forward. Projection. Size[2] = n. Trans. Axial*n. Axial; CL_SAFE_CALL( cl. Enqueue. NDRange. Kernel(commands, forward. Projection. Kernel, 3, NULL, forward. Projection. Size, NULL, 0, NULL) ); cl. Finish(commands); //… }

Előrevetítés: problématér felosztása __kernel void forward. Projection( const int n. Axial, const int n. Trans. Axial, const int resolution, const float detector. Area, float 4 volume. Min, float 4 volume. Max, __global struct module* modules, __global float* measured. Lors, __global float* estimated. Lors, __global float* volume. Buffer ){ int i. Pair = get_global_id(0); int i. Axial 1 = get_global_id(1) % n. Axial; int i. Trans. Axial 1 = get_global_id(1) / n. Axial; int i. Axial 2 = get_global_id(2) % n. Axial; int i. Trans. Axial 2 = get_global_id(2) / n. Axial; int lor. Index = get. Lor. Index( i. Pair, i. Axial 1, i. Axial 2, i. Trans. Axial 1, i. Trans. Axial 2, n. Axial, n. Trans. Axial // … y_estimated kiszámítása if ( y_estimated != 0. 0 f ) estimated. Lors[lor. Index] = measured. Lors[lor. Index ] / y_estimated; else estimated. Lors[lor. Index ] = 0. 0 f; } );

Előrevetítés: problématér felosztása // kitalálhatunk más indexelést is… int get. Lor. Index( int pair, int axial 1, int axial 2, int trans. Axial 1, int trans. Axial 2, int n. Axial, int n. Trans. Axial ) { int lor. Number = n. Axial*n. Trans. Axial; int u = axial 1 * n. Trans. Axial + trans. Axial 1; int v = axial 2 * n. Trans. Axial + trans. Axial 2; return lor. Number*pair + u*n. Axial*n. Trans. Axial + v; }

LOR-képek nézegetése… CL_SAFE_CALL( cl. Enqueue. Read. Buffer(commands, estimated. Lor. Buffer. GPU, CL_TRUE, 0, sizeof(float) * n. Pairs*n. Axial*n. Trans. Axial, estimated. Lor. Buffer. CPU, 0, NULL)); save. TGA(iteration, n. Axial*n. Trans. Axial, n. Trans. Axial*n. Pairs, estimated. Lor. Buffer. CPU); Az 5 D leképzése 2 D-be, általában nem túl átlátható, de néha segíthet. Helyette: a modulokkal párhuzamos LOR-ok kiválogatása

Panelekre merőleges LOR-ok A LOR-ok közül kiválogatjuk azokat, amelyek végpontjai egymás „tükörképei” - Axiális koordináta megegyezik - Tranzaxiálist tükrözzük a 2 - A számításokhoz nincs rá o 2 szükség, de pl. hibajavításhoz t 2 hasznos t 0 a 0 o 0 n 0

Panelekre merőleges LOR-ok cl_mem parallel. Projection. Buffer. GPU; float parallel. Projection. Buffer. CPU[n. Pairs*n. Axial*n. Trans. Axial]; // Panelpáronként annyi merőleges LOR van, ahány detektor egy panelen void init. Simulation(){ // … parallel. Projection. Buffer. GPU = cl. Create. Buffer(context, CL_MEM_READ_WRITE, sizeof(float) * n. Pairs*n. Axial*n. Trans. Axial, NULL); } void simulation. Step(){ //… CL_SAFE_CALL( cl. Set. Kernel. Arg(parallel. Projection. Kernel, 0, sizeof(int), &n. Axial) ); CL_SAFE_CALL( cl. Set. Kernel. Arg(parallel. Projection. Kernel, 1, sizeof(int), &n. Trans. Axial) ); CL_SAFE_CALL( cl. Set. Kernel. Arg(parallel. Projection. Kernel, 2, sizeof(cl_mem), &estimated. Lor. Buffer. GPU) ); CL_SAFE_CALL( cl. Set. Kernel. Arg(parallel. Projection. Kernel, 3, sizeof(cl_mem), &parallel. Projection. Buffer. GPU) ); size_t parallel. Projection. Size[3]; parallel. Projection. Size[0] = n. Pairs; parallel. Projection. Size[1] = n. Axial; parallel. Projection. Size[2] = n. Trans. Axial; CL_SAFE_CALL( cl. Enqueue. NDRange. Kernel(commands, parallel. Projection. Kernel, 3, NULL, parallel. Projection. Size, NULL, 0, NULL) ); } __kernel void parallel. Projection( const int n. Axial, const int n. Trans. Axial, __global float* input. Lors, __global float* output. Lors ){ int i. Pair = get_global_id(0); int i. Axial = get_global_id(1); int i. Trans. Axial = get_global_id(2); int lor. Index = get. Lor. Index( i. Pair, i. Axial, i. Trans. Axial, n. Trans. Axial-i. Trans. Axial-1, n. Axial, n. Trans. Axial ); int parallel. Index = i. Axial + i. Trans. Axial*n. Axial + i. Pair*n. Axial*n. Trans. Axial; output. Lors[parallel. Index] = input. Lors[lor. Index]; }

Panelekre merőleges LOR-ok void simulation. Step(){ //… // parallel projection után CL_SAFE_CALL( cl. Enqueue. Read. Buffer(commands, parallel. Projection. Buffer. GPU, CL_TRUE, 0, sizeof(float) * n. Pairs*n. Axial*n. Trans. Axial, parallel. Projection. Buffer. CPU, 0, NULL)); save. TGA(iteration, n. Axial, n. Trans. Axial*n. Pairs, parallel. Projection. Buffer. CPU); //… }

Honnan szerezzük a mérést? 1. Veszünk/építünk egy PET-készüléket, és betoljuk a laborba. . . 2. Feltételezzük, hogy jó a modellünk, és a saját előrevetítővel állítjuk elő a „mérést” – A valóságban is valami hasonló történik, az igazi mérés a várható érték (amit mi is közelítünk)+némi zaj

Honnan szerezzük a mérést? • A kiindulási volume-ot arra állítjuk, amit szeretnénk rekonstruálni – A valóságban ezt nyilván nem ismerjük // Pl. gömb (órán ezt használjuk): for ( unsigned int x = 0; x < volume. Size[0]; ++x ) for ( unsigned int y = 0; y < volume. Size[1]; ++y ) for ( unsigned int z = 0; z < volume. Size[2]; ++z ) { if ( (x-resolution/2)*(x-resolution/2)+(y-resolution/2)*(y-resolution/2)+(z-resolution/2)*(zresolution/2) < 3*(resolution/5) ) volume. Data[ x + volume. Size[1]*y + volume. Size[0]*volume. Size[1]*z ] = 1. 0 f; else volume. Data[ x + volume. Size[1]*y + volume. Size[0]*volume. Size[1]*z ] = 0. 0 f; } // vagy fájlból betöltés: void load. Volume(const char* file. Name){ FILE* data. File = fopen(file. Name, "rb"); char* magic. Num = new char[2]; fread(magic. Num, sizeof(char), 2, data. File); if('V' == magic. Num[0] && 'F' == magic. Num[1]){ int volume. Size[3]; fread(volume. Size, sizeof(int), 3, data. File); volume. Data = new float[volume. Size[0] * volume. Size[1] * volume. Size[2]]; fread(volume. Data, sizeof(float), volume. Size[0] * volume. Size[1] * volume. Size[2], data. File); resolution = volume. Size[0]; // most köbös volume-ot feltételezünk } else { std: : cout << "Can't open volume file" << file. Name << std: : endl; } }

„Mérésgenerátor” mód __kernel void forward. Projection( // … int measurement. Generation ){ // … if ( measurement. Generation ){ measured. Lors[lor. Index] = y_estimated; return; } if ( 0. 0 f != y_estimated ) estimated. Lors[lor. Index] = measured. Lors[lor. Index] / y_estimated; else estimated. Lors[lor. Index] = 0. 0 f; }

Mérés generálás //… // Volume felöltése a „páciens” adataival (betöltés, procedurális generálás stb. ) //… // előrevetítjük, „mérésgenerátor” módban measurement. Generation = 1; CL_SAFE_CALL( cl. Set. Kernel. Arg(forward. Projection. Kernel, 0, sizeof(int), &n. Axial) ); //… CL_SAFE_CALL( cl. Set. Kernel. Arg(forward. Projection. Kernel, 10, sizeof(cl_mem), &measurement. Generation) ); size_t forward. Projection. Size[3]; forward. Projection. Size[0] = n. Pairs; forward. Projection. Size[1] = forward. Projection. Size[2] = n. Trans. Axial*n. Axial; CL_SAFE_CALL( cl. Enqueue. NDRange. Kernel(commands, forward. Projection. Kernel, 3, NULL, forward. Projection. Size, NULL, 0, NULL) ); cl. Finish(commands); // konstansra állítjuk a kezdeti volume-ot, ebből indul majd a rekonstrukció („normál” módban) measurement. Generation = 0; volume. Size[0] = volume. Size[1] = volume. Size[2] = resolution; volume. Data = new float[resolution * resolution]; for ( unsigned int x = 0; x < resolution; ++x ) for ( unsigned int y = 0; y < resolution; ++y ) for ( unsigned int z = 0; z < resolution; ++z ) { volume. Data[ x +resolution*y +resolution*z ] = 1. 0 f; } CL_SAFE_CALL( cl. Enqueue. Write. Buffer(commands, volume. Buffer. GPU, CL_TRUE, 0, sizeof(float) * resolution, volume. Data, 0, NULL)); • Ezzel véget ért a keretprogram

Előrevetítés: volume-metszés __kernel void forward. Projection( /*…*/ ) { // … indexek számítása float y_estimated = 0. 0 f; int 2 i. Module = get. Modules(i. Pair); float 4 z 1 = modules[i. Module. x]. origin + modules[i. Module. x]. trans. Axial*(i. Trans. Axial 1+0. 5 f)/n. Trans. Axial + modules[i. Module. x]. axial*(i. Axial 1+0. 5 f)/n. Axial; float 4 z 2 = modules[i. Module. y]. origin + modules[i. Module. y]. trans. Axial*(i. Trans. Axial 2+0. 5 f)/n. Trans. Axial + modules[i. Module. y]. axial*(i. Axial 2+0. 5 f)/n. Axial; float 4 dir = z 2 -z 1; float tnear, tfar; if ( intersect. Box( z 1, dir, volume. Min, volume. Max, &tnear, &tfar ) ) { // Ray-marching, y_estimated kiszámítása } // eredmény kiírása }

Melyik LOR metszi a volume-ot? __kernel void forward. Projection( /*…*/ ) { // … indexek számítása float y_estimated = 0. 0 f; int 2 i. Module = get. Modules(i. Pair); float 4 z 1 = modules[i. Module. x]. origin + modules[i. Module. x]. trans. Axial*(i. Trans. Axial 1+0. 5 f)/n. Trans. Axial + modules[i. Module. x]. axial*(i. Axial 1+0. 5 f)/n. Axial; float 4 z 2 = modules[i. Module. y]. origin + modules[i. Module. y]. trans. Axial*(i. Trans. Axial 2+0. 5 f)/n. Trans. Axial + modules[i. Module. y]. axial*(i. Axial 2+0. 5 f)/n. Axial; float 4 dir = z 2 -z 1; float tnear, tfar; if ( intersect. Box( z 1, dir, volume. Min, volume. Max, &tnear, &tfar ) ) { // piros = kísérletezgetés, rekonstrukciós kódnak nem része estimated. Lors[lor. Index] = 1. 0 f; return; } else { estimated. Lors[lor. Index] = 0. 0 f; retun; } }

A volume AABB vetülete a panelekre • A normálissal párhuzamos LOR-ok: Panelpárok • Próbáljuk ki más AABB-re!

Előrevetítés: y_estimated float get. Intensity(const float 4 p, const int resolution, __global float* intensity. Buffer){ int x = p. x * resolution; int y = p. y * resolution; int z = p. z * resolution; if(x > resolution -1 || x < 0) return(0. 0 f); if(y > resolution -1 || y < 0) return(0. 0 f); if(z > resolution -1 || z < 0) return(0. 0 f); return intensity. Buffer[x + y * resolution + z * resolution]; } __kernel void forward. Projection( /*…*/ ) { // … if ( intersect. Box( z 1, dir, volume. Min, volume. Max, &tnear, &tfar ) ) { float G = -detector. Area*detector. Area * dot(modules[i. Module. x]. n, dir)*dot(modules[i. Module. y]. n, dir) / (2. 0 f*M_PI*dot(dir, dir)); float 4 start = z 1+tnear*dir; float 4 end = z 1+tfar*dir; float 4 step = (end - start) / resolution; float dl = length(step); float 4 voxel = start; for ( int i = 0; i < resolution; ++i ) { float x = get. Intensity( (voxel - volume. Min) / (volume. Max - volume. Min), resolution, volume. Buffer ); y_estimated += G*x*dl; voxel += step; } } // eredmény kiírása }

A volume vetülete a detektorokra • measurement. Generation mód be, a merőleges vetítést vizsgáljuk __kernel void forward. Projection( /*…*/ ) { // … if ( intersect. Box( z 1, dir, volume. Min, volume. Max, &tnear, &tfar ) ) { float G = -detector. Area*detector. Area * dot(modules[i. Module. x]. n, dir)*dot(modules[i. Module. y]. n, dir) / (2. 0 f*M_PI*dot(dir, dir)); float 4 start = z 1+tnear*dir; float 4 end = z 1+tfar*dir; float 4 step = (end - start) / resolution; float dl = length(step); float 4 voxel = start; for ( int i = 0; i < resolution; ++i ) { float x = get. Intensity( (voxel - volume. Min) / (volume. Max - volume. Min), resolution, volume. Buffer ); // geom. faktorral nem szorzunk, a vonalintegrált jelenítjük meg // dl sem érdekes, a II-os LOR-oknál úgyis megegyezik y_estimated += x; voxel += step; } } // eredmény kiírása }

A volume vetülete a detektorokra Vonalintegrál • Mi történik, ha tologatjuk, átméretezzük az AABB-t?

Önálló feladat • Szimuláljuk a bináris tomográfiát • Oda írjunk 1 -et, ahol a LOR-nak van metszete az objektummal Vonalintegrál Bináris

Visszavetítés: problématér felosztása // Host: void simulation. Step(){ //… CL_SAFE_CALL( CL_SAFE_CALL( CL_SAFE_CALL( cl. Set. Kernel. Arg(back. Projection. Kernel, cl. Set. Kernel. Arg(back. Projection. Kernel, 0, 1, 2, 3, 4, 5, 6, 7, 8, sizeof(int), &n. Pairs) ); sizeof(int), &n. Axial) ); sizeof(int), &n. Trans. Axial) ); sizeof(int), &resolution) ); sizeof(cl_float 4), &volume. Min) ); sizeof(cl_float 4), &volume. Max) ); sizeof(cl_mem), &module. Buffer. GPU) ); sizeof(cl_mem), &estimated. Lor. Buffer. GPU) ); sizeof(cl_mem), &volume. Buffer. GPU) ); size_t back. Projection. Size[3]; back. Projection. Size[0] = back. Projection. Size[1] = back. Projection. Size[2] = resolution; CL_SAFE_CALL( cl. Enqueue. NDRange. Kernel(commands, back. Projection. Kernel, 3, NULL, back. Projection. Size, NULL, 0, NULL) ); cl. Finish(commands); }

Visszavetítés: problématér felosztása __kernel void back. Projection( const int n. Pairs, const int n. Axial, const int n. Trans. Axial, const int resolution, float 4 volume. Min, float 4 volume. Max, __global struct module* modules, __global float* estimated. Lors, __global float* volume. Buffer ){ int 4 i. Voxel = (int 4)(get_global_id(0), get_global_id(1), get_global_id(2), 0); int linear. Index = i. Voxel. x + i. Voxel. y*resolution + i. Voxel. z*resolution; // … voxel érték kiszámítása if ( denominator > 0. 0 f ) { volume. Buffer[linear. Index] *= numerator / denominator; } else // ha nem lát rá egy LOR sem a voxelre { volume. Buffer[linear. Index] = 0. 0 f; } }

Visszavetítés: voxel a térben __kernel void back. Projection( const int n. Pairs, const int n. Axial, const int n. Trans. Axial, const int resolution, float 4 volume. Min, float 4 volume. Max, __global struct module* modules, __global float* estimated. Lors, __global float* volume. Buffer ){ //… float 4 voxel = (float 4)(i. Voxel. x, i. Voxel. y, i. Voxel. z, 0); voxel /= resolution; // [0, 1] voxel = (volume. Max-volume. Min)*voxel + volume. Min; // [volume. Min, Volume. Max] //… }

Végiggyalogolunk a távolabbi detektoron __kernel void back. Projection( /*…*/ ){ //… float numerator = 0. 0 f; float denominator = 0. 0 f; for ( int i. Pair = 0; i. Pair < n. Pairs; ++i. Pair ) { int 2 i. Module = get. Modules(i. Pair); struct module primary = modules[i. Module. x]; // far detector struct module secondary = modules[i. Module. y]; // near detector bool switch. Detectors = false; if ( dot(voxel-primary. origin, primary. n) < dot(voxel-secondary. origin, secondary. n) ) { switch. Detectors = true; struct module tmp; tmp = primary; primary = secondary; secondary = tmp; } float s. Trans. Axial = n. Trans. Axial / dot(secondary. trans. Axial, secondary. trans. Axial); float s. Axial = n. Axial / dot(secondary. axial, secondary. axial); for ( int p = 0; p < n. Axial; ++p ) for ( int q = 0; q < n. Trans. Axial; ++q ) { //… }

Vetület, LOR hozzájárulása for ( int p = 0; p < n. Axial; ++p ) for ( int q = 0; q < n. Trans. Axial; ++q ) { float 4 z 1 = primary. origin + primary. axial*(p+0. 5 f)/n. Axial + primary. trans. Axial*(q+0. 5 f)/n. Trans. Axial; float 4 dv = voxel - z 1; // using similar triangle for projection float t = dot(secondary. n, secondary. origin-z 1) / dot(secondary. n, dv); float 4 z 2 = z 1 + dv*t; float fr = dot(z 2 -secondary. origin, secondary. axial) * s. Axial; float fs = dot(z 2 -secondary. origin, secondary. trans. Axial) * s. Trans. Axial; if ( 0 <= fr && fr < n. Axial && 0 <= fs && fs < n. Trans. Axial ) { int r = (int)fr; int s = (int)fs; int lor. Index; if ( switch. Detectors ) { lor. Index = get. Lor. Index( i. Pair, r, p, s, q, n. Axial, n. Trans. Axial ); } else { lor. Index = get. Lor. Index( i. Pair, p, r, q, s, n. Axial, n. Trans. Axial ); } float y_div_yw = estimated. Lors[lor. Index]; float ddv = length(dv); float Alv = dot(primary. n, dv) / (ddv*ddv); numerator += y_div_yw * Alv; denominator += Alv; } } }

Voxelek meghatározottsága: térszög • • A voxelekből a teljes detektorgyűrű mekkora térszöget fed le A voxelből induló fotonpárok detektálási valószínűsége __kernel void back. Projection( /*…*/ ){ //… int lor. Cnt = 0; for ( int i. Pair = 0; i. Pair < n. Pairs; ++i. Pair ) { //… for ( int p = 0; p < n. Axial; ++p ) for ( int q = 0; q < n. Trans. Axial; ++q ) { //… // mekkora térszögben látszik a LOR a voxelből denominator += Alv; } //… // a térszögek összegét írjuk ki volume. Buffer[linear. Index] = denominator; }

Voxelek meghatározottsága: térszög YZ XZ • Tologassuk, méretezzük át az AABB-t! Mi történik? XY

Önálló feladat • A térszög helyett a voxelen átmenő LOR-ok számát nézzük meg • Van-e különbség az előzőhöz képest?

Rekonstrukció • Írjunk vissza minden kísérleti változást (szedjük ki a piros részeket)

Önálló/szorgalmi feladat • Változtassuk a volume AABB-t, a kiindulási volume-hoz ne nyúljunk – Változik-e a rekonstrukció eredménye? – Mit kell ahhoz csinálnunk, hogy változzon? • A referencia LOR-halmaz előállítása során megadott térfogat lógjon ki a rekonstrukció által használt AABB-ből – Mi történik? Miért?