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

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] }