Cerebral Blood Flow Measurement (Unarranged)

/*****************************************************************
* program: TC_CBF_std.c
*
* purpose: generate CBF standard file
*
* author: Ta-Cheng Chang
*
* modified: July 11, 2008
*
****************************************************************
/#include "cbook.h"

int TC_CBF_std(void)

{
FILE *fin;
char temp[64];
char * directory = (char*)malloc(512);
char * tempdic1 = (char*)malloc(512);
char * tempdic2 = (char*)malloc(512);
char * tempdic3 = (char*)malloc(512);
char * tempstring = (char*)malloc(512);
char * fname = (char*)malloc(512);
char * stddic = (char*)malloc(64);
int tp, i, j, k, choice, indx;
float lambda, Kvalue, m, W;
float * Flowtable = (float*)malloc(4 * 256);
m = 1; // Net effect diffusion limitation (0~1); default = 1
W = 100; // Tissue mass; default = 100 [g]

clr;
move (12,5); bold;
printf ("Instruction:\n"); move (12,6); bold; printf ("STEP1 "); norm; printf ("Enter animal file directory to import date"); move (19,7); printf ("Folder contains Time.txt, DPM.txt, Intensity.txt\n"); move (12,8); bold; printf ("STEP2 "); norm; printf ("Choose radioisotope concentration standard\n"); move (12,9); bold; printf ("STEP3 "); norm; printf ("Intensity vs. blood flow file (ff.std) is generated at the same directory\n"); move (19,10); printf("Blood flow unit: ml/min/100g"); move (12,11); printf ("\n");
move (12,12); bold; printf ("Enter animal file directory: ex: "); norm; printf ("/usr/users/IMAGES/CBF-DATA/0304164"); move (12,14);
/***************************************************** * Import Data * *****************************************************/ /************** Import Time File ********************/ INQUIRY: scanf ("%s", &directory);
strcpy (tempdic1, &directory); strcpy (tempstring, "/Time.txt"); fname = strcat (tempdic1, tempstring); fin = fopen (fname, "r"); strcpy (tempdic1, &directory);
//fin = fopen ("/usr/users/IMAGES/CBF-DATA/0304164/T0304164.txt", "r"); if ( fin == 0 ) { move (12,13); printf ("No time file! Please re-enter directory!\n"); move (12,14); goto INQUIRY; } fscanf (fin, "%s", temp); tp = atoi (temp); temp[0] = 0;
fscanf (fin, "%s", temp); lambda = atof (temp); temp[0] = 0;
float * tdatam = (float*)malloc(4 * tp); float * tdatas = (float*)malloc(4 * tp); for (i = 0; i < tp; i++) { fscanf (fin, "%s" ,temp); tdatam[i] = atof (temp); //Minute tdatas[i] = atof (temp) * 60; //Second } fclose (fin); temp[0] = 0;
/***************** Import DPM File ***********************/ strcpy (tempdic2, tempdic1); strcpy (tempstring, "/DPM.txt"); fname = strcat (tempdic1, tempstring); fin = fopen (fname, "r"); //fin = fopen ("/usr/users/IMAGES/CBF-DATA/0304164/DPM0304164.txt", "r"); if ( fin == 0 ) { move (12,13); printf ("No DPM file! Please re-enter directory!\n"); move (12,14); goto INQUIRY; }
float * dpmdata = (float*)malloc(4 * tp); float * Ca = (float*)malloc(4 * tp); for (i = 0; i < tp; i++) { fscanf (fin, "%s", temp); dpmdata[i] = atof (temp); Ca[i] = atof (temp) / 2220 / 0.02; //Measured unit:[DPM/l];wanted unit:[nCi/ml].Sample size:20ml.1[nCi]=2220[DPM]. } fclose (fin); temp[0] = 0;
/****************** Import Standard ***********************/ move (12,16); bold; printf ("Choose standard file:\n"); move (12,17); bold; printf ("1) "); norm; printf ("set1\n"); move (12,18); bold; printf ("2) "); norm; printf ("set2\n"); move (12,19); bold; printf ("3) "); norm; printf ("146A\n"); move (12,20); bold; printf ("4) "); norm; printf ("146E\n"); printf ("\n"); move (12,22); bold; printf ("Enter Choice: "); norm;
INQUIRY2: move (26,22);
scanf ("%d", &choice); switch ( choice ) { case 1: stddic = "set1.txt"; break; case 2: stddic = "set2.txt"; break; case 3: stddic = "146A.txt"; break; case 4: stddic = "146E.txt"; break; default: goto INQUIRY2; break; } strcpy(tempstring, "/usr/users/IMAGES/STD/"); fname = strcat (tempstring, stddic); fin = fopen (fname, "r"); if ( fin == 0 ) { move (12,23); printf ("No standard file! Please check!\n"); goto INQUIRY2; }
i = 0; float * tempstd = (float*)malloc(4 * 30); while(fscanf (fin, "%s", temp) != EOF ) { tempstd[i] = atof (temp); i++; } int stdpt = i; float * stdata = (float*)malloc(4 * stdpt); for (i = 0; i < stdpt; i++) { stdata[i] = tempstd[i]; } free(tempstd); fclose (fin); temp[0] = 0;
/**************** Import Standard Intenstity *****************/ strcpy (tempdic3, tempdic2); strcpy (tempstring, "/Intensity.txt"); fname = strcat (tempdic2, tempstring); fin = fopen (fname, "r");
//fin = fopen ("/usr/users/IMAGES/CBF-DATA/0304164/I0304164.txt", "r"); if ( fin == 0 ) { move (12,13); printf ("No intensity file! Please re-enter directory!\n"); move (12,14); goto INQUIRY; } int stdint[stdpt]; for (i = 0; i < stdpt; i++) { fscanf (fin, "%s", temp); stdint[i] = atoi (temp); } fclose (fin); temp[0] = 0;
/******************* Printing ************************/ move (12,24); printf ("Total number of time points is %d\n", tp); move (12,25); printf ("Lambda is %f [ml/g]\n", lambda); move (12,26); printf("\n"); move (12,27); bold; printf("Sample Time[sec] Ca[nCi/ml] Ca[DPM]\n"); norm; for (i = 0; i < tp; i++) { move (12,28 + i); printf(" %d %f %f %f\n", i + 1, tdatas[i], Ca[i], dpmdata[i]); } indx = 28 + i; move (12,indx + 1); bold; printf("Count Intensity Ci[nCi/g]\n"); norm; for (i = 0; i < stdpt; i++) { move (12,indx + 2 + i); printf(" %d %d %f\n", i + 1, stdint[i], stdata[i]); } indx = indx + 2 + i; printf("\n");
/************************************************************ * Standard Intensity vs Standard Concentration Curve Fitting * *************************************************************/ float intervall; float * cfconint = (float*)malloc(4 * 256);
/******************* 0 to the 1st Point *********************/ intervall = stdata[0] / stdint[0]; for (i = 0; i <= stdint[0] - 1; i++) { cfconint[i] = intervall * i; }
/****** From the 1st Point to Rest Points *******************/ for (i = 1; i <= stdpt - 1; i++) { cfconint[stdint[i - 1]] = stdata[i - 1]; if (stdint[i - 1] == stdint[i]) { stdint[i] = stdint[i] + 1; } intervall = (stdata[i] - stdata[i - 1]) / (stdint[i] - stdint[i - 1]); for (j = 1; j <= (stdint[i] - stdint[i - 1]); j++) { cfconint[stdint[i - 1] + j] = intervall * j + stdata[i - 1]; } }
/********** From the last Point to 255 ***********************/ for (i = 1; i <= 255-stdint[stdpt - 1]; i++) { cfconint[stdint[stdpt - 1] + i] = intervall * i + stdata[stdpt - 1]; }
/********************** Printing curve fitting result ***************************** bold; printf("Intensity Concentration[nCi/g]\n"); norm; for (i = 0; i <= 255; i++) { printf(" %d %f\n", i, cfconint[i]); } printf("\n"); printf("\n");*/
/************************************************************ * Generating Blood Flow Table * *************************************************************/ /****** Generate K table (0~0.35 in 10^-5 precision) ********/ int Ctlength, index; float tempk, temp2; float * tempCt = (float*)malloc(4 * 35100);
for (i = 0; i <= 35000; i++) { tempk = i * 0.00001; temp2 = 0; for (j = 1; j <= tp - 1; j++) { temp2 = temp2 + lambda * tempk * Ca[j] * exp(-tempk * (45 - tdatas[j])) * (tdatas[j] - tdatas[j - 1]); } tempCt[i] = temp2; if (i > 0) { if (tempCt[i] < tempCt[i - 1]) { break; } } } Ctlength = i; float * Ct = (float*)malloc(4 * (Ctlength - 1)); //Resize Ct for (j = 0; j <= Ctlength - 1; j++) { Ct[j] = tempCt[j]; } free(tempCt);
/****** Approach K value in 10^-6 precision; generate bloodflow table ********/ for (i = 0; i <= 255; i++) { for (j = 0; j <= Ctlength -1; j++) { if (cfconint[i] < index =" j;" temp2 =" 0;" j =" 0;" kvalue =" (0.00001)" k =" 1;" temp2 =" temp2">= cfconint[i]) { Flowtable[i] = Kvalue * 60 * lambda * W / m; // Bloodflow unit: ml/min/100g; default m = 1; W = 100 [g] break; } } }
/****** Generate blood flow file at the same directory **********/ strcpy (tempstring, "/ff.std"); fname = strcat (tempdic3, tempstring); fin = fopen (fname, "w"); for (i = 0; i <= 254; i++) { fprintf(fin, "%d %f\n", i, Flowtable[i]); } fprintf(fin, "0 0\n"); fclose (fin);
/*****************************************************************/ free(tempdic1); free(tempdic2); free(tempdic3); free(tempstring); free(tdatam); free(tdatas); free(dpmdata); free(Ca); free(stdata); free(Flowtable); move (12,indx + 1); printf("Type 'm' then press enter to back to main menu: "); move (61,indx + 1); scanf ("%d",&k); }

Region Growing

To be continued....

LU Decomposition

LU decomposition is a method to solve matrix questions. By finding the "pivot", the original matrix can be transfer to an upper triangular matrix. The unknown matrix then can be easily solved by simple algebra calculation. Take the following matrix for example, the c3Z=C, b2Y+b3Z=B, a1X+a2Y+a3Z=A. Put Z=C/c3 into above equations, the X and Y can be easily solved.

Tcl/Tk Code:
## LU Decomposition. Find X if AX=B by LU decomposition. After calculation X=B.
for {set i 0} {$i <= 1} {incr i 1} {
for {set k [expr {$i + 1}]} {$k <= 2} {incr k 1} {
set temp [expr {double($A($k,$i)) / double($A($i,$i))}]
for {set l [expr {$i + 1}]} {$l <= 2} {incr l 1} {
set A($k,$l) [expr {$A($k,$l) - ($A($i,$l) * $temp)}]
}
set B($k) [expr {$B($k) - (($B($i)) * $temp)}]
}
}
set B(2) [expr {$B(2) / $A(2,2)}]
for {set i 1} {$i >= 0} {incr i -1} {
for {set k [expr {$i + 1}]} {$k <= 2} {incr k 1} {
set B($i) [expr {$B($i) - $A($i,$k) * $B($k)}]
}
set B($i) [expr {$B($i) / $A($i,$i)}]
}

Best Fit Circle Using Least Squares Method

For best fitting a circle with a set of data points, the least squares method has been considered as one of the most useful method. The least squares method is based on minimizing the sum of mean square distance from the circle perimeter to data points. By given n points (xi, yi), 1<=i<=n, the objective function can be defined as

where di is the distance from the point (xi, yi) to the perimeter. The equation of circle is

where (a, b) is the circle center and R is the radius. Then the di is

The equation of the least squares method can be rewritten as

where B=-2a, C=-2b, and D=a^2+b^2-R^2. Next, differentiating F with respect to B, C, D yields a system of linear equations

B, C, and D can be solved by solving matrix

Let

The equation of matrix calculation can be rewritten as

The X can be solved by calculating the equation

where A^-1 is the inverse matrix of A. The matrix X can be solved using LU decomposition to get B, C, and D. The x coordinate of the center (a) equals to -B/2, the y coordinate of the center (b) equals to -C/2, and the radius (R) equals to

Reference:

N. Chernov and C. Lesort, Least squares fitting of circles, J. Math. Imaging & Vision 23, 2005, 239-252.

Code:

Tcl/Tk

## Circular regression (least squares method) by given the coordinates of a set of data points.

## Return the center points (XC , YC), and radius R. => CircleFit X Y returns XC YC R

proc CircleFit { X Y } {

set n [llength $X]

catch {unset XX ; unset YY ; unset XY ; unset XXY ; unset YYY ; unset XXX ; unset XYY}
set sumX 0 ; set sumY 0 ; set sumXX 0 ; set sumYY 0 ; set sumXY 0
set sumXXY 0 ; set sumYYY 0 ; set sumXXX 0 ; set sumXYY 0
set sumXXplusYY 0 ; set sumXXYplusYYY 0 ; set sumXXXplusXYY 0
for {set i 0} {$i < $n} {incr i 1} { set XX [lappend XX [expr {[lindex $X $i] * [lindex $X $i]}]] set YY [lappend YY [expr {[lindex $Y $i] * [lindex $Y $i]}]] set XY [lappend XY [expr {[lindex $X $i] * [lindex $Y $i]}]] set XXY [lappend XXY [expr {[lindex $X $i] * [lindex $X $i] * [lindex $Y $i]}]] set YYY [lappend YYY [expr {[lindex $Y $i] * [lindex $Y $i] * [lindex $Y $i]}]] set XXX [lappend XXX [expr {[lindex $X $i] * [lindex $X $i] * [lindex $X $i]}]] set XYY [lappend XYY [expr {[lindex $X $i] * [lindex $Y $i] * [lindex $Y $i]}]] set sumX [expr {$sumX + [lindex $X $i]}] set sumY [expr {$sumY + [lindex $Y $i]}] set sumXX [expr {$sumXX + [lindex $XX $i]}] set sumYY [expr {$sumYY + [lindex $YY $i]}] set sumXY [expr {$sumXY + [lindex $XY $i]}] set sumXXY [expr {$sumXXY + [lindex $XXY $i]}] set sumYYY [expr {$sumYYY + [lindex $YYY $i]}] set sumXXX [expr {$sumXXX + [lindex $XXX $i]}] set sumXYY [expr {$sumXYY + [lindex $XYY $i]}] } set nsumXXplusYY -[expr {$sumXX + $sumYY}] set nsumXXYplusYYY -[expr {$sumXXY + $sumYYY}] set nsumXXXplusXYY -[expr {$sumXXX + $sumXYY}] set A(0,0) $sumX ; set A(0,1) $sumY ; set A(0,2) $n set A(1,0) $sumXY ; set A(1,1) $sumYY ; set A(1,2) $sumY set A(2,0) $sumXX ; set A(2,1) $sumXY ; set A(2,2) $sumX set B(0) $nsumXXplusYY ; set B(1) $nsumXXYplusYYY ; set B(2) $nsumXXXplusXYY ## Gaussian elimination. Find X if AX=B by LU decomposition. After calculation X=B. for {set i 0} {$i <= 1} {incr i 1} { for {set k [expr {$i + 1}]} {$k <= 2} {incr k 1} { set temp [expr {double($A($k,$i)) / double($A($i,$i))}] for {set l [expr {$i + 1}]} {$l <= 2} {incr l 1} { set A($k,$l) [expr {$A($k,$l) - ($A($i,$l) * $temp)}] } set B($k) [expr {$B($k) - (($B($i)) * $temp)}] } } set B(2) [expr {$B(2) / $A(2,2)}] for {set i 1} {$i >= 0} {incr i -1} {
for {set k [expr {$i + 1}]} {$k <= 2} {incr k 1} { set B($i) [expr {$B($i) - $A($i,$k) * $B($k)}] } set B($i) [expr {$B($i) / $A($i,$i)}] } ## Calculate the center and radius of the circle set XC [expr {-(0.5 * $B(0))}] set YC [expr {-(0.5 * $B(1))}] set R [expr {sqrt((pow($B(0),2) + pow($B(1),2)) / 4 - $B(2))}] return [list $XC $YC $R] }

A New Face of this Blog

Finally, I decided to switch this blog space to my experimental "nerd base". The first reason is that I think it can be a good review to track all the works I've done. The other reason is I hope it can help people who may meet the same problems and is looking for solutions. I'll only post articles related to my research, study, and work. Please feel free to leave a message if you have any idea or better solution, or, just make a comment^^. Enjoy~~~~