// this file contains: mri.cpp, mri.h, and brainmri128.h /* =============================== cut here ============================= */ /* ============================== mri.cpp ============================ */ /* =============================== cut here ============================= */ /******************************************************************* mri.cpp Revision 1.0 5/4/2002 under construction ********************************************************************/ /* Copyright (C) 2002 Doug Neubauer * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include #include #include #include #include #include #include "mri.h" #include "brainmri128.h" // const double DTIMEROTATING = (2.0e-7); // const double DTIMEROTATINGREADOUT = (2.010e-7); const double DTIMEROTATINGREADOUT = (SAMPLEPERIOD / 20.0); const double DTIMEROTATING = (1.0e-6); const double DTIMEROTATINGFID = (4.0e-6); #if 1 commandList Brains::pulseSequence[] = { {LABORATORYFRAME, 1.5, 0.0, 0.0 }, /* Tesla */ {EQUILIBRIUMMAGNETIZATION, 0.0, 0.0, 0.0 }, {ROTATINGFRAME,DTIMEROTATING, 0.0, 0.0 }, /* ns, speed up simulation */ {RESETMAGNETICFIELD,0.0, 0.0, 0.0 }, /* speed up simulation (force tRecovery) */ {RESETMRISTATE,-GRADIENTPERVOXELY, 0.0, 0.0 }, /* Y phase, cumulativeDelay = -GRADIENTPERVOXELY */ {DOLOOP,NUMBEROFDISPLAYVOXELSY, 0.0, 0.0 }, // {DOLOOP,1, 0.0, 0.0 }, {RESETMAGNETICFIELD,0.0, 0.0, 0.0 }, /* speed up simulation (force tRecovery) */ // {ENABLEPRINT,1.0e-6, 0.0, 0.0 }, {SETGRADIENT, 0.0, 0.0, GRADIENTPERVOXELZ }, {BLOCHEQUATION, SLICEDURATION, 0.0, B1MAGNITUDE }, /* rf90 */ // {DUMPVOXELS,0.0, 0.0, 0.0 }, {SETGRADIENTY , GRADIENTPERVOXELX, dGRADIENTPERVOXELY, 0.0 }, {BLOCHEQUATION, PHASEDURATION, 0.0, 0.0 }, /* phase */ {SETGRADIENT,0.0, 0.0, 0.0 }, /* B0 only, no gradients */ {ROTATINGFRAME,DTIMEROTATINGFID, 0.0, 0.0 }, /* ns, speed up simulation */ {BLOCHEQUATION, ECHOTIME / 2.0 - SLICEDURATION - PHASEDURATION, 0.0, 0.0 }, {ROTATINGFRAME,DTIMEROTATING, 0.0, 0.0 }, /* ns, speed up simulation */ {SETGRADIENT, 0.0, 0.0, GRADIENTPERVOXELZ }, {BLOCHEQUATION, SLICEDURATION * 1.0, 0.0, PHASEADJUSTMENTVALUE * 2.0 * B1MAGNITUDE }, /* rf180 */ {SETGRADIENT,0.0, 0.0, 0.0 }, /* B0 only, no gradients */ {ROTATINGFRAME,DTIMEROTATINGFID, 0.0, 0.0 }, /* ns, speed up simulation */ {BLOCHEQUATION, ECHOTIME / 2.0 - SLICEDURATION / 2.0 - READOUTTIME / 2.0, 0.0, 0.0 }, {ROTATINGFRAME,DTIMEROTATING, 0.0, 0.0 }, /* ns, speed up simulation */ {SETGRADIENT ,GRADIENTPERVOXELX, 0.0, 0.0 }, {ROTATINGFRAME,DTIMEROTATINGREADOUT, 0.0, 0.0 }, /* ns, speed up simulation */ {BLOCHEQUATION, READOUTTIME + DTIMEROTATINGREADOUT, NUMBEROFDISPLAYVOXELSX, 0.0 }, /* readout */ {ROTATINGFRAME,DTIMEROTATING, 0.0, 0.0 }, /* ns, speed up simulation */ // {SETGRADIENT,0.0, 0.0, 0.0 }, /* B0 only, no gradients */ // {BLOCHEQUATION, 2.0e-3, 0.0, 0.0 }, // {DISABLEPRINT,0.0, 0.0, 0.0 }, // {DUMPVOXELS,0.0, 0.0, 0.0 }, {ENDLOOP,0.0, 0.0, 0.0 }, }; #endif #if 0 commandList Brains::pulseSequence[] = { {LABORATORYFRAME, 1.5, 0.0, 0.0 }, /* Tesla */ {EQUILIBRIUMMAGNETIZATION, 0.0, 0.0, 0.0 }, {RESETMAGNETICFIELD,0.0, 0.0, 0.0 }, /* speed up simulation (force tRecovery) */ {RESETMRISTATE,-GRADIENTPERVOXELY, 0.0, 0.0 }, /* Y phase, cumulativeDelay = -GRADIENTPERVOXELY */ {DOLOOP,1, 0.0, 0.0 }, {RESETMAGNETICFIELD,0.0, 0.0, 0.0 }, /* speed up simulation (force tRecovery) */ {ENABLEPRINT,1.0e-7, 0.0, 0.0 }, {SETGRADIENT, 0.0, 0.0, GRADIENTPERVOXELZ }, {BLOCHEQUATION, SLICEDURATION, 0.0, B1MAGNITUDE }, /* rf90 */ // {DUMPVOXELS,0.0, 0.0, 0.0 }, {SETGRADIENT,0.0, 0.0, 0.0 }, /* B0 only, no gradients */ {BLOCHEQUATION, 2.0e-4, 0.0, 0.0 }, {DISABLEPRINT,0.0, 0.0, 0.0 }, // {DUMPVOXELS,0.0, 0.0, 0.0 }, {ENDLOOP,0.0, 0.0, 0.0 }, }; #endif int main( ) { Brains brainTop; // test(); srand(0); brainTop.RunBrain(); return(0) ; } #if 0 char *initializeVoxelsArray[NUMBEROFVOXELSY * 1] = { //01234567 x //000111122223333 "................", // y 15 z = 0 ".....kkkkkk.....", // y 14 z = 0 "...kkmWWWWmkk...", // y 13 z = 0 "..kmmWWWWWWmmk..", // y 12 z = 0 "..kmWWGGGGWWmk..", // y 11 z = 0 "..kWWGGGGGGWWk..", // y 10 z = 0 "..kWWGGGGGGWWk..", // y 9 z = 0 "..kWWWGWWGWWWk..", // y 8 z = 0 "..kWWcGWWGcWWk..", // y 7 z = 0 "..kWWcGWWGcWWk..", // y 6 z = 0 "..kWWccGGccWWk..", // y 5 z = 0 "..kWWWccccWWWk..", // y 4 z = 0 "..kWWWWWWWWWWk..", // y 3 z = 0 "...kkWWWWWWkk...", // y 2 z = 0 ".....kkkkkk.....", // y 1 z = 0 "................", // y 0 z = 0 }; #endif #if 0 char *initializeVoxelsArray[NUMBEROFVOXELSY * 1] = { //01234567 x //000111122223333 "wwwwFwwwwwwwwwww", // y 7 z = 0 "wwwwwwwwwwwwwwww", // y 6 z = 0 "wwwwwwwwwwwwwwww", // y 5 z = 0 "wwwwwwwwwwwwwwww", // y 4 z = 0 "wwwwwwwwwwwwwwww", // y 3 z = 0 "wwwwwwwwwwwwwwww", // y 2 z = 0 "wwwwFwwwwwwwwwww", // y 1 z = 0 "wwwwFwwwwwwwwwww", // y 0 z = 0 }; #endif #if 0 char *initializeVoxelsArray[NUMBEROFVOXELSY * 4] = { //01234567 x "...W", // y 3 z = 0 "....", // y 2 "....", // y 1 "....", // y 0 "..F.", // y 3 z = 2 "F..F", // y 2 "F.F.", // y 1 ".F..", // y 0 ".F.F", // y 3 z = 1 ".F.F", // y 2 ".F.F", // y 1 ".F.F", // y 0 "....", // y 3 z = 3 "FFFF", // y 2 "....", // y 1 "FFFF", // y 0 }; #endif Brains::Brains () { int i, j,jj, k,ii; char *s ; char ch ; double spinValue ; double t2Value; fprintf(stderr,"DTIMEROTATINGREADOUT:%g, DTIMEROTATING:%g\n",DTIMEROTATINGREADOUT,DTIMEROTATING); hBar = 6.6260754e-34 / (2.0 * PIEVALUE); /* joules-second */ boltzmannConstant = 1.380657e-23; /* joules/K */ temperature = 37.0 + 273.15; /* 310 Kelvin (body temperature) */ printf("temperature: %g Kelvin\n",temperature); printf("boltzmannConstant: %g joules/K\n",boltzmannConstant); printf("hBar: %g joules-second\n",hBar); printFlag = 0; printTime = 0.0 ; dprintTime = 1.0e-6; sampleCount = 0; sampleTime = 0.0 ; currentTime = 0.0; localFrequency = 0.0 ; dLocalFrequency = 0.0; rfCoilSensitivity = 0.0; laboratoryFlag = 1; brain = new Voxel **[NUMBEROFVOXELSZ] ; for (i = 0; i < NUMBEROFVOXELSZ; i++) { brain[i] = new Voxel *[NUMBEROFVOXELSY] ; for (j = 0; j < NUMBEROFVOXELSY; j++) { brain[i][j] = new Voxel[NUMBEROFVOXELSX]; } /* jForLoop */ } /* iForLoop */ for (k = 0; k < NUMBEROFVOXELSZ; k++) { for (j = 0; j < NUMBEROFVOXELSY; j++) { jj = j / (NUMBEROFVOXELSY / NUMBEROFDISPLAYVOXELSY); // s = initializeVoxelsArray[(NUMBEROFDISPLAYVOXELSY - jj - 1) + ((k & 3) * 4)] ; s = initializeVoxelsArray[(NUMBEROFVOXELSY - j - 1) + ((k & 3) * 4)] ; for (i = 0; i < NUMBEROFVOXELSX; i++) { ii = i / (NUMBEROFVOXELSX / NUMBEROFDISPLAYVOXELSX); // ch = s[ii]; ch = s[i]; switch (ch) { case '.': /* background, empty */ // printf("empty[z%d][y%d][x%d]\n",k,j,i); spinValue = SPINSPERVOXEL * 0.0 ; t2Value = T2EMPTY; break; case 'c': /* csf */ spinValue = SPINSPERVOXEL * 1.0 ; t2Value = 0.2 / SPEEDUPFACTOR; break; case 'G': /* gray matter */ spinValue = SPINSPERVOXEL * 0.84 ; t2Value = T2GRAY; break; case 'W': /* white matter */ spinValue = SPINSPERVOXEL * 0.72 ; t2Value = T2WHITE; break; case 'F': /* fat */ spinValue = SPINSPERVOXEL * 0.9 ; t2Value = T2FAT; // printf("fat[z%d][y%d][x%d]\n",k,j,i); break; case 'm': /* muscle/skin */ spinValue = SPINSPERVOXEL * 0.9 ; /* ? */ t2Value = 0.040 / SPEEDUPFACTOR; break; case 's': /* skin */ spinValue = SPINSPERVOXEL * 0.9 ; /* ? */ t2Value = 0.050 / SPEEDUPFACTOR; /* ? */ break; case 'k': /* skull */ spinValue = SPINSPERVOXEL * 0.9 ; /* ? */ t2Value = 0.020 / SPEEDUPFACTOR; /* ? */ break; case 'g': /* glial matter */ spinValue = SPINSPERVOXEL * 0.9 ; /* ? */ t2Value = 0.080 / SPEEDUPFACTOR; /* ? */ break; case 't': /* connective */ spinValue = SPINSPERVOXEL * 0.9 ; /* ? */ t2Value = 0.030 / SPEEDUPFACTOR; /* ? */ break; case 'w': /* water */ // printf("water[z%d][y%d][x%d]\n",k,j,i); spinValue = SPINSPERVOXEL * 1.0 ; t2Value = T2WATER; break; default: fprintf(stderr,"error unknown %c %d\n",ch,ch); exit(1); break; } /* switch */ brain[k][j][i].InitializeVoxel(spinValue,t2Value,GYROMAGNETICVALUE); } /* iForLoop */ } /* jForLoop */ } /* kForLoop */ } void Brains::RunBrain () { static const int NUMBEROFCOMMANDS = (sizeof(Brains::pulseSequence) / sizeof(commandList)); static const NUMBEROFEMBEDDEDLOOPS = (10); double dTime ; int i, j, k; double larmorFrequency ; double magneticField0Z; double currentGradientY = 0.0; int command; int commandIndex; int loopStack[NUMBEROFEMBEDDEDLOOPS] ; int loopIndex = 0; int loopCount; rfCoilSensitivity = 0.000234885; /* Tesla/amp, C=B1/A, assume: 1kv, 1kw, 1kohm = 1amp */ currentTime = 0.0; sampleTime = 0.0; sampleCount = 0; laboratoryFlag = 1; printTime = 0.0 ; dprintTime = 0.0; StatusDisplay(); // fprintf(stderr,"NUMBEROFCOMMANDS: %d\n",NUMBEROFCOMMANDS); for (commandIndex = 0; commandIndex < NUMBEROFCOMMANDS; commandIndex++) { command = pulseSequence[commandIndex].command ; fprintf(stderr,"commandIndex:%d\n",commandIndex); switch (command) { case LABORATORYFRAME: dTime = 1.0e-9; /* 1ns, expected by low pass filter design */ // dTime = 2.0e-9; /* 1ns, expected by low pass filter design */ printf("LABORATORYFRAME: dTime: %g\n",dTime); magneticField0Z = pulseSequence[commandIndex].value[0]; /* Tesla */ larmorFrequency = GYROMAGNETICVALUE * magneticField0Z; dLocalFrequency = larmorFrequency * dTime ; // for local oscillators printf("magneticField0[z]: %g Tesla\n",magneticField0Z); printf("larmorFrequency: %g Hertz\n", larmorFrequency/ (2.0 * PIEVALUE) ); laboratoryFlag = 1; break; case ROTATINGFRAME: dTime = pulseSequence[commandIndex].value[0]; printf("ROTATINGFRAME: dTime: %g\n",dTime); magneticField0Z = 0.0; /* Tesla */ larmorFrequency = 0.0; dLocalFrequency = 0.0; localFrequency = 0.0 ; laboratoryFlag = 0; break; case ENABLEPRINT: printFlag = 1; dprintTime = pulseSequence[commandIndex].value[0]; /* seconds */ break; case DISABLEPRINT: printFlag = 0; break; case DOLOOP: loopIndex++; if (loopIndex >= NUMBEROFEMBEDDEDLOOPS) { fprintf(stderr,"loopStack overflow, %d\n",commandIndex); exit(1); } /* if */ loopStack[loopIndex] = commandIndex ; loopCount = (int )pulseSequence[commandIndex].value[0]; break; case ENDLOOP: printf("ENDLOOP: %d\n",loopCount); fprintf(stderr,"ENDLOOP: %d\n",loopCount); if (loopIndex <= 0) { fprintf(stderr,"loopStack under flow, %d\n",commandIndex); exit(1); } /* if */ loopCount-- ; if (loopCount <= 0) { loopIndex-- ; break; } /* if */ commandIndex = loopStack[loopIndex]; break; case RESETMRISTATE: currentGradientY = pulseSequence[commandIndex].value[0]; break; case EQUILIBRIUMMAGNETIZATION: for (k = 0; k < NUMBEROFVOXELSZ; k++) { for (j = 0; j < NUMBEROFVOXELSY; j++) { for (i = 0; i < NUMBEROFVOXELSX; i++) { brain[k][j][i].SetEquilibriumMagnetization(magneticField0Z, hBar,boltzmannConstant,temperature); } /* iForLoop */ } /* jForLoop */ } /* kForLoop */ break; case DUMPVOXELS: for (k = 0; k < NUMBEROFVOXELSZ; k++) { for (j = 0; j < NUMBEROFVOXELSY; j++) { for (i = 0; i < NUMBEROFVOXELSX; i++) { printf("Voxel[z%d][y%d][x%d]: ",k,j,i); brain[k][j][i].DumpVoxel(); } /* iForLoop */ } /* jForLoop */ } /* kForLoop */ break; case RESETMAGNETICFIELD: for (k = 0; k < NUMBEROFVOXELSZ; k++) { for (j = 0; j < NUMBEROFVOXELSY; j++) { for (i = 0; i < NUMBEROFVOXELSX; i++) { brain[k][j][i].ResetBulkMagnetization(); /* hack: to speed up simulation */ } /* iForLoop */ } /* jForLoop */ } /* kForLoop */ break; case SETGRADIENT: SetMagneticField0(magneticField0Z,pulseSequence[commandIndex].value[0],pulseSequence[commandIndex].value[1],pulseSequence[commandIndex].value[2]); break; case SETGRADIENTY: SetMagneticField0(magneticField0Z,pulseSequence[commandIndex].value[0],currentGradientY,pulseSequence[commandIndex].value[2]); printf("SETGRADIENTY: currentGradientY:%g\n",currentGradientY); currentGradientY += pulseSequence[commandIndex].value[1]; break; case BLOCHEQUATION: sampleCount = (int )pulseSequence[commandIndex].value[1]; /* number of samples */ sampleTime = currentTime + SAMPLEPERIOD; /* synchronize */ ExecuteRf( pulseSequence[commandIndex].value[0], dTime, pulseSequence[commandIndex].value[2]); if (sampleCount != 0) { fprintf(stderr,"sampleCount under flow, %d, currentTime:%g\n",sampleCount,currentTime); } /* if */ break; default: fprintf(stderr,"illegal: %d\n",command); break; } /* switch */ } /* commandIndexForLoop */ } void Voxel::InitializeVoxel(double numberOfSpins1, double t2Initialize ,double gyroMagneticRatio1) { numberOfSpins = numberOfSpins1; gyroMagneticRatio = gyroMagneticRatio1; t1 = 0.26 ; t2 = t2Initialize; spinDensity = 0.0; magneticFieldB0Z = 0.0; bulkMagnetizationX = 0.0; bulkMagnetizationY = 0.0; bulkMagnetizationZ = 0.0; equilibriumMagnetizationZ = 0.0; } double Voxel::fieldVariations[16] = { +5.0e-6 * VARIATIONSMULTIPLIER, +4.28e-6 * VARIATIONSMULTIPLIER, +3.57e-6 * VARIATIONSMULTIPLIER, +2.86e-6 * VARIATIONSMULTIPLIER, +2.14e-6 * VARIATIONSMULTIPLIER, +1.43e-6 * VARIATIONSMULTIPLIER, +0.714e-6 * VARIATIONSMULTIPLIER, +0.0e-6 * VARIATIONSMULTIPLIER, -5.0e-6 * VARIATIONSMULTIPLIER, -4.28e-6 * VARIATIONSMULTIPLIER, -3.57e-6 * VARIATIONSMULTIPLIER, -2.86e-6 * VARIATIONSMULTIPLIER, -2.14e-6 * VARIATIONSMULTIPLIER, -1.43e-6 * VARIATIONSMULTIPLIER, -0.714e-6 * VARIATIONSMULTIPLIER, -0.0e-6 * VARIATIONSMULTIPLIER, }; void Voxel::SetEquilibriumMagnetization (double magneticFieldValue, double hBar,double boltzmannConstant,double temperature) { double magneticMoment, dEnergy, spinRatio ; magneticVariation = magneticFieldValue * fieldVariations[rand() & 0x0f]; magneticMoment = 0.5 * gyroMagneticRatio * hBar ; dEnergy = gyroMagneticRatio * hBar * magneticFieldValue ; spinRatio = exp( dEnergy / (boltzmannConstant * temperature) ) ; spinDensity = numberOfSpins * (spinRatio - 1.0) / (spinRatio + 1.0) ; // printf("spinDensity: %g per %s spinRatio:%20.10g\n",spinDensity,VOXELUNITS,spinRatio); equilibriumMagnetizationZ = spinDensity * magneticMoment; // printf("larmorFrequency: %g Hertz\n", gyroMagneticRatio * magneticFieldValue/ (2.0 * PIEVALUE)); // bulkMagnetizationX = 1.0e-12; // bulkMagnetizationY = 2.0e-12; // bulkMagnetizationZ = equilibriumMagnetizationZ; /* assumed at equilibrium */ } void Voxel::ResetBulkMagnetization () { bulkMagnetizationX = 0.0; bulkMagnetizationY = 0.0; bulkMagnetizationZ = equilibriumMagnetizationZ; /* hack: force to equilibrium */ // printf("ResetBulkMagnetization: %g\n",equilibriumMagnetizationZ); } void Voxel::PrintVoxel () { #if 0 printf("numberOfSpins: %g per %s\n",numberOfSpins,VOXELUNITS); printf("gyroMagneticRatio: %g radians/(second * Tesla)\n",gyroMagneticRatio); printf("spinDensity: %g per %s\n",spinDensity,VOXELUNITS); printf("equilibriumMagnetizationZ: %g amps/meter\n",equilibriumMagnetizationZ * VOXELSIZE2CUBICMETERS); #endif printf("Mx: %g\n",bulkMagnetizationX); printf("My: %g\n",bulkMagnetizationY); printf("Mz: %g\n",bulkMagnetizationZ); // printf("b0: %g\n",magneticFieldB0Z); } void Voxel::DumpVoxel () { printf("Mx:%g My:%g Mz:%g (MzEq:%g)\n",bulkMagnetizationX,bulkMagnetizationY,bulkMagnetizationZ,equilibriumMagnetizationZ); } void test () { int i,j,ch; double x,y; char *s; for (y = 0.0,j = 0; j < 217; j++, y += 0.59447) { if (y < 1.0) { continue; } /* if */ y -= 1.0; s = initializeVoxelsArray[j]; for (x = 0.0,i = 0; i < 181; i++,x += 0.712707) { if (x < 1.0) { continue; } /* if */ x -= 1.0; ch = s[i] ; printf("%c",ch); } /* iForLoop */ printf("\n"); } /* jForLoop */ #if 0 char convertTable[] = ".cGWFmskgt"; int i,j,k,ch; FILE *filePointer ; filePointer = fopen("e:\\temp\\brainweb1.exe","rb") ; for (k = 0; k < 181; k++) { // printf("z: %d\n",k); for (j = 0; j < 217; j++) { for (i = 0; i < 181; i++) { ch = fgetc(filePointer) ; if (ch == EOF) { fprintf(stderr,"end of file\n"); exit(1); } /* if */ if (ch > 9) { fprintf(stderr,"ch: %d\n",ch); ch = 0; } /* if */ ch = convertTable[ch] ; if (k == 90) { printf("%c",ch); } /* if */ } /* iForLoop */ if (k == 90) { printf("\n"); } /* if */ } /* jForLoop */ // printf("\n"); } /* forLoop */ double k ; double dFrequency,frequency,gz; for (k = -10; k < 10; k++) { gz = 10.0e-6 * 10.0 * k; frequency = gz * GYROMAGNETICVALUE / (2.0 * PIEVALUE) ; dFrequency = GYROMAGNETICVALUE * 10.0e-6 * 5.0 / (2.0 * PIEVALUE); printf("k:%g\n",k); printf(" %g\n",gz); printf(" %g\n",frequency - dFrequency); printf(" %g\n",frequency); printf(" %g\n\n",frequency + dFrequency); } /* kForLoop */ printf("Lamar frequency: %g\n",1.5 * GYROMAGNETICVALUE / (2.0 * PIEVALUE)); #endif exit(0); } void Brains::StatusDisplay () { char *title = "Status" ; char *titleHeading = "status" ; FILE *filePointer; int i,j,k; double dx,dy,x1; char ch; filePointer = fopen("status.html","w") ; fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer," \n"); fprintf(filePointer," \n"); fprintf(filePointer," \n"); fprintf(filePointer," %s\n",title); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"
\n"); fprintf(filePointer,"

\n"); fprintf(filePointer,"%s

\n",titleHeading); fprintf(filePointer,"\n"); fprintf(filePointer,"pie:%g
\n",PIEVALUE); fprintf(filePointer,"FIELDOFVIEW: x:%g y:%g z:%g (mm)
\n",FIELDOFVIEWX,FIELDOFVIEWY,FIELDOFVIEWZ); fprintf(filePointer,"NUMBEROFVOXELS: x:%d y:%d z:%d
\n",NUMBEROFVOXELSX,NUMBEROFVOXELSY,NUMBEROFVOXELSZ); fprintf(filePointer,"NUMBEROFDISPLAYVOXELS: x:%d y:%d z:%d
\n",NUMBEROFDISPLAYVOXELSX,NUMBEROFDISPLAYVOXELSY,NUMBEROFDISPLAYVOXELSZ); fprintf(filePointer,"gyroMagneticRatio: %g (radians/(second * Tesla), hydrogen) %g (Hz/Tesla)
\n",GYROMAGNETICVALUE,GYROMAGNETICVALUE / (2.0 * PIEVALUE)); fprintf(filePointer,"VOXELSIZE: %g %s
\n",VOXELSIZE, VOXELUNITS); fprintf(filePointer,"SPINSPERVOXEL: %g
\n",SPINSPERVOXEL); fprintf(filePointer,"T2EMPTY: %g
\n",T2EMPTY); fprintf(filePointer,"T2FAT: %g
\n",T2FAT); fprintf(filePointer,"T2WATER: %g
\n",T2WATER); fprintf(filePointer,"ECHOTIME: %g
\n",ECHOTIME); fprintf(filePointer,"GRADIENTXVALUE: %g (Tesla/mm)
\n",GRADIENTXVALUE); fprintf(filePointer,"GRADIENTPERVOXELX: %g (Tesla/Voxel)
\n",GRADIENTPERVOXELX); fprintf(filePointer,"BANDWIDTH: %g (Hz)
\n",BANDWIDTH); fprintf(filePointer,"SAMPLEPERIOD:%g (seconds), sampleFrequency:%g (Hz)
\n",SAMPLEPERIOD,BANDWIDTH); fprintf(filePointer,"READOUTTIME: %g (seconds)
\n",READOUTTIME); fprintf(filePointer,"PHASEDURATION: %g (seconds)
\n",PHASEDURATION); fprintf(filePointer,"GRADIENTYMAXIMUM: %g (Tesla/mm)
\n",GRADIENTYMAXIMUM); fprintf(filePointer,"GRADIENTPERVOXELY: %g (Tesla/Voxel, maximum)
\n",GRADIENTPERVOXELY); fprintf(filePointer,"dGRADIENTPERVOXELY: %g (Tesla/Voxel, minimum)
\n",dGRADIENTPERVOXELY); fprintf(filePointer,"GRADIENTZVALUE: %g (Tesla/mm)
\n",GRADIENTZVALUE); fprintf(filePointer,"GRADIENTPERVOXELZ: %g (Tesla/Voxel)
\n",GRADIENTPERVOXELZ); fprintf(filePointer,"SLICETHICKNESS: %g (mm)
\n",SLICETHICKNESS); fprintf(filePointer,"NUMBEROFLOBES: %g
\n",NUMBEROFLOBES); fprintf(filePointer,"EFFECTIVELENGTH: %g (seconds)
\n",EFFECTIVELENGTH); fprintf(filePointer,"SLICEDURATION: %g (seconds)
\n",SLICEDURATION); fprintf(filePointer,"B1MAGNITUDE: %g (Tesla), %g (Hz)
\n",B1MAGNITUDE,B1MAGNITUDE * GYROMAGNETICVALUE / (2.0 * PIEVALUE)); fprintf(filePointer,"Hz/Voxel: x: %g == %g
\n",1.0 / (SAMPLEPERIOD * NUMBEROFVOXELSX), GRADIENTPERVOXELX * GYROMAGNETICVALUE / (2.0 * PIEVALUE)); fprintf(filePointer,"Hz/Voxel: Y: %g Hz (assume sampleRateY = 1Hz) (check: %g == %g)
\n",1.0 / ( NUMBEROFVOXELSY), PIEVALUE, PHASEDURATION * GRADIENTYMAXIMUM * GYROMAGNETICVALUE * FIELDOFVIEWY / (double )NUMBEROFDISPLAYVOXELSY); fprintf(filePointer,"

\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); for (i = 0; i < NUMBEROFVOXELSX; i++) { fprintf(filePointer,"\n",i); } /* iForLoop */ fprintf(filePointer,"\n"); for (j = NUMBEROFVOXELSY-1; j >= 0 ; j--) { dy = (GRADIENTPERVOXELY * PHASEDURATION * GYROMAGNETICVALUE / ((double )NUMBEROFDISPLAYVOXELSY * PIEVALUE) ) * (double )(j - VOXELCENTERY); fprintf(filePointer,"\n"); fprintf(filePointer,"\n", j); for (i = 0; i < NUMBEROFVOXELSX; i++) { dx = GRADIENTPERVOXELX * GYROMAGNETICVALUE / (2.0 * PIEVALUE) * (double )(i - VOXELCENTERX); fprintf(filePointer,"\n",dy ,dx); } /* iForLoop */ } /* jForLoop */ fprintf(filePointer,"
K Space (numbers in Hertz)
  X
Y  		%d
%d%.3f
%.0f
\n"); fprintf(filePointer,"

\n"); fprintf(filePointer,"

\n");
    for (k = 0; k < NUMBEROFVOXELSZ; k++) {
        fprintf(filePointer,"\n    z:%d\n",k);
        for (i = 0; i < NUMBEROFVOXELSX; i++) {
            fprintf(filePointer,"%d",i%10);
        } /* iForLoop */
        fprintf(filePointer,"\n");
        for (j = NUMBEROFVOXELSY-1; j  >= 0 ; j--) {
            for (i = 0; i < NUMBEROFVOXELSX; i++) {
                x1 = brain[k][j][i].ReadT2();
                if (x1 == T2FAT) {
                    ch = 'f';
                } /* if */
                else if (x1 == T2WATER) {
                    ch = 'w';
                } /* elseIf */     
                else {
                    ch = '.';
                } /* else */
                fprintf(filePointer,"%c",ch);
            } /* iForLoop */
            fprintf(filePointer,"  // y:%d\n",j);
        } /* jForLoop */
    } /* kForLoop */
    fprintf(filePointer,"
\n"); fprintf(filePointer,"

\n"); fprintf(filePointer,"\n"); fprintf(filePointer,"\n"); fclose(filePointer); } void Voxel::ExecuteBloch (double dTime ,double magneticField1X,double magneticField1Y, double &signalX,double &signalY) { // adapted from: WWW.eecs.umich.edu keywords: dnoll blochsim // alternate form where rotations are explicitly calculated and carried out // on the magnetization vector - this is a more stable solution // alternate (from Runge) method is faster and less sensitive to dt ... // (at least 100x faster than Eulers.. dn) // // Compute sines & cosines of field angles: // Theta = angle w.r.t positive z axis // Phi = angle w.r.t positive x axis // Psi = angle w.r.t transformed positive x axis // double x1,x2 ; double magneticFieldMagnitude ; double transverseMagnitude ; double cosineTheta ; double sineTheta; double cphi,sphi,spsi,cpsi; double Mx0,My0,Mz0; double Mx1,My1,Mz1; double localmagneticField0Z; localmagneticField0Z = magneticFieldB0Z * gyroMagneticRatio * dTime; magneticField1X *= gyroMagneticRatio * dTime; magneticField1Y *= gyroMagneticRatio * dTime; x1 = (magneticField1X * magneticField1X) + (magneticField1Y * magneticField1Y); x2 = (localmagneticField0Z * localmagneticField0Z) + x1; // Magnitude^2 of applied field if (x2 != 0.0) { magneticFieldMagnitude = sqrt( x2 ); // Magnitude of applied field Mx0 = bulkMagnetizationX; My0 = bulkMagnetizationY; Mz0 = bulkMagnetizationZ; transverseMagnitude = sqrt(x1); // Magnitude of transverse applied field cosineTheta = localmagneticField0Z / magneticFieldMagnitude; // cos(theta) sineTheta = sqrt(1.0 - (cosineTheta * cosineTheta)); // sin(theta) > 0 cphi = 1.0; if (transverseMagnitude != 0.0) { cphi = magneticField1X / transverseMagnitude; // cos(phi) } /* if */ sphi = sqrt(1.0 - (cphi * cphi)); // sin(phi) if (magneticField1Y < 0.0) { sphi = -sphi; } /* if */ cpsi = cos(magneticFieldMagnitude); // cos(psi) spsi = sin(magneticFieldMagnitude); // sin(psi) Mx1 = cphi*(cosineTheta*(cpsi*(cosineTheta*(sphi*My0+cphi*Mx0)-sineTheta*Mz0) + spsi*(cphi*My0-sphi*Mx0))+sineTheta*(cosineTheta*Mz0+sineTheta*(sphi*My0+cphi*Mx0))) - sphi*(-spsi*(cosineTheta*(sphi*My0+cphi*Mx0)-sineTheta*Mz0) + cpsi*(cphi*My0-sphi*Mx0)); My1 = sphi*(cosineTheta*(cpsi*(cosineTheta*(sphi*My0+cphi*Mx0)-sineTheta*Mz0) + spsi*(cphi*My0-sphi*Mx0))+sineTheta*(cosineTheta*Mz0+sineTheta*(sphi*My0+cphi*Mx0))) + cphi*(-spsi*(cosineTheta*(sphi*My0+cphi*Mx0)-sineTheta*Mz0) + cpsi*(cphi*My0-sphi*Mx0)); Mz1 = cosineTheta*(cosineTheta*Mz0+sineTheta*(sphi*My0+cphi*Mx0)) - sineTheta*(cpsi*(cosineTheta*(sphi*My0+cphi*Mx0)-sineTheta*Mz0) + spsi*(cphi*My0-sphi*Mx0)); } /* if */ else { Mx1 = bulkMagnetizationX; My1 = bulkMagnetizationY; Mz1 = bulkMagnetizationZ; } /* else */ // relaxation effects: bulkMagnetizationX = Mx1 - dTime * Mx1 / t2; bulkMagnetizationY = My1 - dTime * My1 / t2; bulkMagnetizationZ = Mz1 - dTime * (Mz1 - equilibriumMagnetizationZ) / t1; #if 0 bulkMagnetizationX = Mx1 - (dTime * Mx1 / t2) * (1.0 - 0.5 * dTime/t2); bulkMagnetizationY = My1 - (dTime * My1 / t2) * (1.0 - 0.5 * dTime/t2); bulkMagnetizationZ = Mz1 - dTime * (Mz1 - equilibriumMagnetizationZ) / t1; #endif signalX += bulkMagnetizationX ; signalY += bulkMagnetizationY ; #if 0 localmagneticField0Z *= gyroMagneticRatio * dTime; magneticField1X *= gyroMagneticRatio * dTime; magneticField1Y *= gyroMagneticRatio * dTime; bulkMagnetizationX += bulkMagnetizationY * localmagneticField0Z; // bulkMagnetizationX -= bulkMagnetizationZ * magneticField1Y; bulkMagnetizationY -= bulkMagnetizationX * localmagneticField0Z; bulkMagnetizationY += bulkMagnetizationZ * magneticField1X; // bulkMagnetizationZ += bulkMagnetizationX * magneticField1Y; bulkMagnetizationZ -= bulkMagnetizationY * magneticField1X; bulkMagnetizationX -= dTime * bulkMagnetizationX / t2; bulkMagnetizationY -= dTime * bulkMagnetizationY / t2; bulkMagnetizationZ -= dTime * (bulkMagnetizationZ - equilibriumMagnetizationZ) / t1; localmagneticField0Z *= gyroMagneticRatio * dTime; magneticField1X *= gyroMagneticRatio * dTime; magneticField1Y *= gyroMagneticRatio * dTime; bulkMagnetizationX += bulkMagnetizationY * localmagneticField0Z; bulkMagnetizationX -= bulkMagnetizationZ * magneticField1Y; bulkMagnetizationX -= dTime * bulkMagnetizationX / t2; bulkMagnetizationY += bulkMagnetizationZ * magneticField1X; bulkMagnetizationY -= bulkMagnetizationX * localmagneticField0Z; bulkMagnetizationY -= dTime * bulkMagnetizationY / t2; bulkMagnetizationZ += bulkMagnetizationX * magneticField1Y; bulkMagnetizationZ -= bulkMagnetizationY * magneticField1X; bulkMagnetizationZ -= dTime * (bulkMagnetizationZ - equilibriumMagnetizationZ) / t1; double dbulkMagnetizationX ; double dbulkMagnetizationY ; double dbulkMagnetizationZ ; dbulkMagnetizationX = gyroMagneticRatio * ((bulkMagnetizationY * localmagneticField0Z) - (bulkMagnetizationZ * magneticField1Y)) - bulkMagnetizationX / t2; bulkMagnetizationX += dbulkMagnetizationX * dTime; dbulkMagnetizationY = gyroMagneticRatio * ((bulkMagnetizationZ * magneticField1X) - (bulkMagnetizationX * localmagneticField0Z)) - bulkMagnetizationY / t2; bulkMagnetizationY += dbulkMagnetizationY * dTime; dbulkMagnetizationZ = gyroMagneticRatio * ((bulkMagnetizationX * magneticField1Y) - (bulkMagnetizationY * magneticField1X)) - (bulkMagnetizationZ - equilibriumMagnetizationZ) / t1; bulkMagnetizationZ += dbulkMagnetizationZ * dTime; #endif } void Brains::SetMagneticField0 (double magneticField0Z, double gradientX, double gradientY, double gradientZ ) { int i,j,k; double dx,dy,dz; // printf("magneticField0Z:%g \n",magneticField0Z); for (k = 0; k < NUMBEROFVOXELSZ; k++) { dz = gradientZ * (double )(k - VOXELCENTERZ); for (j = 0; j < NUMBEROFVOXELSY; j++) { dy = gradientY * (double )(j - VOXELCENTERY); for (i = 0; i < NUMBEROFVOXELSX; i++) { dx = gradientX * (double )(i - VOXELCENTERX); brain[k][j][i].SetMagneticField0(magneticField0Z + dz + dy + dx); // printf("Voxel[z%d][y%d][x%d]: %g Tesla (Hz/Voxel: x: %g) (dy:%g)\n",k,j,i,magneticField0Z + dz + dy + dx, dx * GYROMAGNETICVALUE / (2.0 * PIEVALUE),dy); } /* iForLoop */ } /* jForLoop */ } /* kForLoop */ } void Brains::ExecuteRf ( double stopTime, double dTime,double magneticField1) { int i,j,k; static double lastSignalX = 0.0 ; static double lastSignalY = 0.0 ; double signalX = 0.0; double signalY = 0.0; double localOscillatorX = 0.0 ; double localOscillatorY = 0.0 ; double demodulatedSignalX,filteredSignalX; double demodulatedSignalY,filteredSignalY; double magneticField1X,magneticField1Y; int toggle; double bfirX0,bfirX1,bfirX2,bfirX3,bfirX4,bfirX5,bfirX6,bfirX7,bfirX8,bfirX9,bfirX10,bfirX11,bfirX12,bfirX13,bfirX14; double bfirY0,bfirY1,bfirY2,bfirY3,bfirY4,bfirY5,bfirY6,bfirY7,bfirY8,bfirY9,bfirY10,bfirY11,bfirY12,bfirY13,bfirY14; double sincTime,dsincTime,sincPulse; // printf("ExecuteRf stopTime:%g, dTime:%g magneticField1:%g\n",stopTime,dTime,magneticField1); sincTime = -NUMBEROFLOBES * PIEVALUE / 2.0 ; dsincTime = NUMBEROFLOBES * PIEVALUE * dTime / SLICEDURATION; stopTime += currentTime; for ( ; currentTime < stopTime; localFrequency += dLocalFrequency, currentTime += dTime, sincTime += dsincTime) { signalX = 0.0; signalY = 0.0; localOscillatorX = cos(localFrequency) ; localOscillatorY = sin(localFrequency) ; if (sincTime == 0.0) { sincPulse = 1.0; } /* if */ else { sincPulse = sin(sincTime) / sincTime; } /* else */ magneticField1X = magneticField1 * sincPulse * localOscillatorX; magneticField1Y = -magneticField1 * sincPulse * localOscillatorY; // magneticField1Y = 0.0; for (k = 0; k < NUMBEROFVOXELSZ; k++) { for (j = 0; j < NUMBEROFVOXELSY; j++) { for (i = 0; i < NUMBEROFVOXELSX; i++) { brain[k][j][i].ExecuteBloch(dTime , magneticField1X, magneticField1Y,signalX,signalY); } /* iForLoop */ } /* jForLoop */ } /* kForLoop */ /* note: signal = rfCoilSensitivity * dM/dt */ if (laboratoryFlag == 0) { demodulatedSignalX = signalX; demodulatedSignalY = -signalY; filteredSignalX = demodulatedSignalX; filteredSignalY = demodulatedSignalY; } /* if */ else { demodulatedSignalX = ((signalX - lastSignalX) * rfCoilSensitivity / dTime) * localOscillatorX; demodulatedSignalY = -((signalY - lastSignalY) * rfCoilSensitivity / dTime) * localOscillatorX; // is this right? // demodulatedSignalY = -signalX * localOscillatorY; // Why '-' ? Revisit!! lastSignalX = signalX; lastSignalY = signalY; toggle++ ; if ((toggle & 1) != 0) { // num_taps:15 freqL:6.5e+007 sampleFrequency:5e+008 (filter:low pass window:Hamming) filteredSignalX = (-0.00194925 * (bfirX0 + bfirX14))+ (-0.00654288 * (bfirX1 + bfirX13))+ (-0.0130404 * (bfirX2 + bfirX12))+ (-0.00436489 * (bfirX3 + bfirX11))+ (0.0434446 * (bfirX4 + bfirX10))+ (0.13133 * (bfirX5 + bfirX9))+ (0.221468 * (bfirX6 + bfirX8))+ (0.26 * bfirX7); bfirX14 = bfirX13; bfirX13 = bfirX12; bfirX12 = bfirX11; bfirX11 = bfirX10; bfirX10 = bfirX9; bfirX9 = bfirX8; bfirX8 = bfirX7; bfirX7 = bfirX6; bfirX6 = bfirX5; bfirX5 = bfirX4; bfirX4 = bfirX3; bfirX3 = bfirX2; bfirX2 = bfirX1; bfirX1 = bfirX0; bfirX0 = demodulatedSignalX; filteredSignalY = (-0.00194925 * (bfirY0 + bfirY14))+ (-0.00654288 * (bfirY1 + bfirY13))+ (-0.0130404 * (bfirY2 + bfirY12))+ (-0.00436489 * (bfirY3 + bfirY11))+ (0.0434446 * (bfirY4 + bfirY10))+ (0.13133 * (bfirY5 + bfirY9))+ (0.221468 * (bfirY6 + bfirY8))+ (0.26 * bfirY7); bfirY14 = bfirY13; bfirY13 = bfirY12; bfirY12 = bfirY11; bfirY11 = bfirY10; bfirY10 = bfirY9; bfirY9 = bfirY8; bfirY8 = bfirY7; bfirY7 = bfirY6; bfirY6 = bfirY5; bfirY5 = bfirY4; bfirY4 = bfirY3; bfirY3 = bfirY2; bfirY2 = bfirY1; bfirY1 = bfirY0; bfirY0 = demodulatedSignalY; } /* if */ } /* else */ if (currentTime >= printTime) { printTime += dprintTime ; if (printFlag) { // fprintf(stderr,"%g\n",currentTime); printf("{ %g\n",currentTime); printf("sampleCount: %d\n",sampleCount); printf("magneticField1X: %g\n",magneticField1X); printf("magneticField1Y: %g\n",magneticField1Y); printf("sincPulse: %g // dsincTime: %g\n",sincPulse,dsincTime); printf("signalX: %g // volts\n",signalX); printf("demodulatedSignalX: %g // volts\n",demodulatedSignalX); printf("filteredsignalX: %g // volts\n",filteredSignalX); printf("signalY: %g // volts\n",signalY); printf("demodulatedSignalY: %g // volts\n",demodulatedSignalY); printf("filteredsignalY: %g // volts\n",filteredSignalY); brain[VOXELCENTERZ][VOXELCENTERY][VOXELCENTERX].PrintVoxel(); printf("}\n"); } /* if */ } /* if */ if (currentTime >= sampleTime) { sampleTime += SAMPLEPERIOD ; if (sampleCount) { printf("{ %g // sampleTime:%g \n",currentTime,sampleTime); printf("xr: %g // SAMPLE, volts\n",filteredSignalX); printf("xi: %g // SAMPLE, volts\n",filteredSignalY); // printf("xi: %g // SAMPLE, volts\n",0.0); printf("}\n"); sampleCount-- ; if (sampleCount == 0) { break; /* all done with samples */ } /* if */ } /* if */ } /* if */ } /* forLoop */ } /* =============================== cut here ============================= */ /* ============================== end mri.cpp ========================= */ /* =============================== cut here ============================= */ /* =============================== cut here ============================= */ /* ============================== mri.h ============================ */ /* =============================== cut here ============================= */ /******************************************************************* MRI.h Revision 1.0 5/4/2002 under construction ********************************************************************/ /* Copyright (C) 2002 Doug Neubauer * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ /* revisit!: dn 5/6/02: possible artifact Bug at Y = NUMBEROFVOXELSY/2 (i.e. in the middle), due to PHASEADJUSTMENTVALUE ?? */ const double PHASEADJUSTMENTVALUE = (1.11); /* needed? because truncated sinc? */ /* <<<< note: some of the numbers are hack'ed to speed up the simulation >>>> */ const double SPEEDUPFACTOR = (5.0); const double SPEEDUPFACTORZ = (100.0); const double PIEVALUE = (atan(1.0) * 4.0); // const double FIELDOFVIEWX = (250.0); /* mm */ // const double FIELDOFVIEWY = (250.0); /* mm */ // const double FIELDOFVIEWZ = (200.0); /* mm */ const double FIELDOFVIEWX = (128.0); /* mm */ const double FIELDOFVIEWY = (128.0); /* mm */ const double FIELDOFVIEWZ = (1.0); /* mm */ /* set VARIATIONSMULTIPLIER = 2.0 in order to test T2* and de-phasing stuff, otherwise zero */ const double VARIATIONSMULTIPLIER = (0.0); // const double VARIATIONSMULTIPLIER = (2.0); /* note: NUMBEROFDISPLAYVOXELS are "normal" viewable voxels */ const int NUMBEROFDISPLAYVOXELSX = (128); const int NUMBEROFDISPLAYVOXELSY = (128); // const int NUMBEROFDISPLAYVOXELSX = (16); // const int NUMBEROFDISPLAYVOXELSY = (16); // const int NUMBEROFDISPLAYVOXELSX = (4); // const int NUMBEROFDISPLAYVOXELSY = (4); const int NUMBEROFDISPLAYVOXELSZ = (1); // const int NUMBEROFDISPLAYVOXELSZ = (64); /* NUMBEROFVOXELS are "sub-voxels" of NUMBEROFDISPLAYVOXELS for testing T2* and de-phasing stuff, otherwise set equal to NUMBEROFDISPLAYVOXELS */ const int NUMBEROFVOXELSX = (1 * NUMBEROFDISPLAYVOXELSX); const int NUMBEROFVOXELSY = (1 * NUMBEROFDISPLAYVOXELSY); // const int NUMBEROFVOXELSX = (4 * NUMBEROFDISPLAYVOXELSX); // const int NUMBEROFVOXELSY = (4 * NUMBEROFDISPLAYVOXELSY); const int NUMBEROFVOXELSZ = (1 * NUMBEROFDISPLAYVOXELSZ); const int VOXELCENTERZ = (NUMBEROFVOXELSZ / 2); const int VOXELCENTERY = (NUMBEROFVOXELSY / 2); const int VOXELCENTERX = (NUMBEROFVOXELSX / 2); const double VOXELSIZE = (FIELDOFVIEWX * FIELDOFVIEWY * FIELDOFVIEWZ / (NUMBEROFVOXELSX * NUMBEROFVOXELSY * NUMBEROFVOXELSZ)); // cubic millimeters const char *VOXELUNITS = ("mm3"); const double VOXELSIZE2CUBICMETERS = (1.0e9); const double SPINSPERVOXEL = (3.35e19 * 2.0 * VOXELSIZE); // mm3, 2 hydrogens const double GYROMAGNETICVALUE = (2.675e8); // radians/(second * Tesla), hydrogen const double T2EMPTY = (0.0001); const double T2WATER = (2.300 / SPEEDUPFACTOR); const double T2FAT = (0.06 / SPEEDUPFACTOR); const double T2WHITE = (0.069 / SPEEDUPFACTOR); const double T2GRAY = (0.106 / SPEEDUPFACTOR); const double ECHOTIME = (T2GRAY * T2WHITE * log(T2GRAY / T2WHITE) / (T2GRAY - T2WHITE)); const double GRADIENTXVALUE = (10.0e-6); /* Tesla/mm */ const double GRADIENTPERVOXELX = (GRADIENTXVALUE * FIELDOFVIEWX / (double )NUMBEROFVOXELSX); const double BANDWIDTH = (GYROMAGNETICVALUE * GRADIENTXVALUE * FIELDOFVIEWX / (2.0 * PIEVALUE)); const double SAMPLEPERIOD = (1.0 / BANDWIDTH); const double READOUTTIME = (SAMPLEPERIOD * (double )NUMBEROFDISPLAYVOXELSX); const double PHASEDURATION = (READOUTTIME / 2.0); /* seconds */ const double GRADIENTYMAXIMUM = (PIEVALUE * (double )NUMBEROFDISPLAYVOXELSY / (FIELDOFVIEWY * PHASEDURATION * GYROMAGNETICVALUE)); const double GRADIENTPERVOXELY = (GRADIENTYMAXIMUM * FIELDOFVIEWY / (double )NUMBEROFVOXELSY); const double dGRADIENTPERVOXELY = (2.0 * GRADIENTPERVOXELY / (double )NUMBEROFDISPLAYVOXELSY); const double GRADIENTZVALUE = (10.0e-6 * SPEEDUPFACTORZ); /* Tesla/mm */ const double GRADIENTPERVOXELZ = (GRADIENTZVALUE * FIELDOFVIEWZ / (double )NUMBEROFVOXELSZ); const double SLICETHICKNESS = (FIELDOFVIEWZ / (double )NUMBEROFDISPLAYVOXELSZ); const double NUMBEROFSIDELOBES = (1 * 2); /* 1 lobe on each side (speed up simulation) */ // const double NUMBEROFSIDELOBES = (2 * 2); const double NUMBEROFLOBES = (NUMBEROFSIDELOBES + 2.0); const double EFFECTIVELENGTH = (2.0 * PIEVALUE / ( GYROMAGNETICVALUE * SLICETHICKNESS * GRADIENTZVALUE)); const double SLICEDURATION = (NUMBEROFLOBES * EFFECTIVELENGTH); const double B1MAGNITUDE = (1.0 * PIEVALUE / (2.0 * GYROMAGNETICVALUE * EFFECTIVELENGTH)); enum { EQUILIBRIUMMAGNETIZATION, LABORATORYFRAME, ROTATINGFRAME, SETGRADIENT, SETGRADIENTY, DISABLEPRINT, ENABLEPRINT, RESETMAGNETICFIELD, DOLOOP, ENDLOOP, RESETMRISTATE, DUMPVOXELS, BLOCHEQUATION, }; #define NUMBEROFVALUES (3) typedef struct _commandList{ int command ; double value[NUMBEROFVALUES]; } commandList ; class Voxel { public: // accessible by all Voxel(){InitializeVoxel(0.0,0.0,0.0);} void Voxel::InitializeVoxel(double numberOfSpins1, double t2Initialize ,double gyroMagneticRatio1); void Voxel::PrintVoxel (); void Voxel::DumpVoxel (); void Voxel::SetEquilibriumMagnetization (double magneticFieldValue, double hBar,double boltzmannConstant,double temperature); void Voxel::ResetBulkMagnetization (); inline void Voxel::ExecuteBloch (double dTime ,double magneticField1X,double magneticField1Y, double &signalX,double &signalY); void SetMagneticField0(double value){magneticFieldB0Z = value + magneticVariation; } double ReadT2(){return(t2) ;} protected: // accessible by subclass private: // accessible by this class only static double fieldVariations[] ; double t1,t2; double gyroMagneticRatio ; double numberOfSpins; double spinDensity; double magneticFieldB0Z; double magneticVariation; double bulkMagnetizationX ; double bulkMagnetizationY ; double bulkMagnetizationZ ; double equilibriumMagnetizationZ ; }; class Brains { public: // accessible by all void RunBrain (); Brains(); protected: // accessible by subclass private: // accessible by this class only static commandList pulseSequence[]; void Brains::ExecuteRf ( double stopTime, double dTime,double magneticField1); void Brains::SetMagneticField0 (double magneticField0Z, double gradientX, double gradientY, double gradientZ ); void Brains::StatusDisplay (); double temperature ; double boltzmannConstant; double hBar; Voxel ***brain; double printTime ; double dprintTime; int printFlag; double sampleTime ; int sampleCount; double currentTime; double dLocalFrequency,localFrequency; double rfCoilSensitivity; int laboratoryFlag; }; /* =============================== cut here ============================= */ /* ============================== end mri.h ========================= */ /* =============================== cut here ============================= */ /* =============================== cut here ============================= */ /* ============================== brainmri128.h ============================ */ /* =============================== cut here ============================= */ /* brainmri128.h */ char *initializeVoxelsArray[128] = { "................................................................................................................................", "................................................................................................................................", "................................................................................................................................", "................................................................................................................................", "................................................................................................................................", "...............................................sssssssssssssssssssssssssssssss..................................................", "............................................sssssssssssssssmmssssssssssssssssssss...............................................", "........................................sssssssmmmmmmmtmmmtmttmmmmtttttmmmmmssssssssss..........................................", ".....................................sssssmmmttttmsssssssssssssssssssmmmmmtmmttmmmsssssss.......................................", "....................................sssmmmmtttmmssskkkkkkkkkkkkkssssssssssssmmmttmmmmmsssss.....................................", "..................................sssmmtttmmssskkkkkkkkkkkkkkkkkkkkkkkkkkkkksssssmmmtmmmssssss..................................", "................................ssmmmttmmssskkkkkssskskkkssssssssssskkkkkkkkkkkkkssssmtttttmmsss................................", "..............................ssmmmttttmsskkkkkssssskkkkkkkkkssssssssssssskkkkkkkkksssmmttttmmsss...............................", "............................ssmmmttttsskkkkkssscckcccccccGGccckkkkkkkkkkcccckkkkkkkkkkssssmtttmmmss.............................", "...........................smmmtFFmsskkkkkkkkkkccGGGGWWGGGWWGGcsssscccccGGGGGcckkkkkkkkkksssmtttmmmss...........................", "..........................smmmttttsskkkkkkkkkkkcGGGGGWWGGWWWGGccsmmGGGGGGGGGGGGcckkkkkkkkksssmtttmmss...........................", "........................smmmtFFtsskkkkkkkkkccccGGWWWGGWWWWWWWGGGcscGGccWWWGGGGGGGGcckkkkkkkksssmtttmmss.........................", "......................ssmmtFFtmskkkkkkkkccGGGGGcGWWWWGGWWWGGGWWGcccccccGWWWWGGGWWGGGcckkkkkkkksssmtttmms........................", "......................smmttFtmsskkkkkkkcGGGGGGGccGWWWWGWWGGGGGGGcccccGGGGGWWWGGWWWGGcckkkkkkkkksssmttmmms.......................", "....................ssmmtFFtsskkkkkkkcGGGWGGWWGGcGGGWWWWWWWGGGGccccGWWWGGGWWWWWWWGGGGGccckkkkkkkkssttttmms......................", "...................sssmtFttskkkskkkcGGGWWWGGWWGGGcGGWWWWWWWGGGGGccGWWWWWWWWWWWWWGGGGGGGGGGckkkkkkkksstFtmmss....................", ".................sssmtFFtmskkskkkcGGGWWWWWGGWWWGGGcGWWWWWWWGGGGGGcGGGWWWWWWWWWWGcGWGWWWWWGGGGckkskkssmttttsss...................", ".................ssmmtFFmsskkkkkccGGGGWWWWGWWWWWGGGGWWWWWWWWGGGGGcGGGGWWWWWWWWWccGGGWWWWGGGGGGkkskkkssmtFtmss...................", "................ssmtFFtmskkkkkccGGGGGGGWWWWWWWWWWGGGWWWWWWWWWWWGGcWWWWWWWWWWWWWGGGGGGGWWGGGWGGGckkkkkssstFtmsss.................", "...............ssmtFFFsskkkkkcGGGGGccGGWWWWWWWWWWWWWWWWWWGGGGGGGGGWWWWWWWWWWWWWWWWGGGGWWWWWWWWGGckkkkkkssmFtmss.................", "..............sssmtFFtskkkkkkcGGWWGGcGGWWWWWWWWWWWWWWWWWWGGGGGGGcGWWWWWWWWWWWWWWWWWGGGWWWWWWWWGGckkkkkkksmFttmss................", "..............ssmtFtmskkkskkccGGWWGGGGGGWWWWWWWWWWWWWWWWWWWGGGccccGGGGWWWWWWWWWWWWWGGGGWWWWWWWGGGcckkkkkkstFtmmss...............", ".............ssmmFFmsskkkkkcccGGGWWWGGGGWWWWWWWWWWWWWWWWWWWWWWGGcGGGGWWWWWWWWWWWWWWGGGGGGGWWWWGGGGGckkkkkkstFtmmss..............", ".............ssmtFFsskkkkkkccGGGGWWWGGGGWWWWWWWWWWWWWWWWWWGWWWGGcGWWWWWWWWWWWWWWWWWGGGGGGGWWWWWGGGGckkkkkksmttmmss..............", "............ssmtFFtskkskkkcccGGGGGWWWGGWWWWWWWWWWWWWWWWWWGGcGGGGGGWWWWWWWWWWWWWWWWWGGGGGGWWWWWWWWGGGckkkkksssFtmmms.............", "............sstFFtmkkkkkccGGGGGGGGWWWWWWWWWWWWWWWWWWWWGGWWWGGGGGGGGGGGGGWWWWWWWWWWWGGGWWWWWWWWWWGGGGGckkkkkssmFtmmss............", "............sstFFmskkkkkcGGGGGGGGGGWWWWWWWWWWWWWWWWWWWGGGWWWWGGGGGGGGGGWWWWWWWWWWWWWGGWWWWWWWWWGGGGGGGckkkkksstttmms............", "...........ssmFFFsskskkcGGGGGGGGGGGGWWWWWWWWWWWWWWWWWWGGGGGGGGGGccGGGWWWWWWGGWWWWWWWWWWWWWWWWWWWGGGWWGGGkkkkksmtFtsm............", "..........ssmtFFmskskkcGGGWWWWGGGGGWWWWWWWWWWWWWWWWWWWGGGGGGGGGGccGGGWGGGGGGGGWWWWWWWWWWWWWWWWWWWWGGGWGGcckkksstFFmsss..........", "..........smmtFtmkkskkGGGWWWWWGGGGWWWWWWWWWWWWWWWWWWWWWGGGGGGGGccGGGGGGGGGGGGGWWWWWWWWWWWWWWWWWWGGGGGGGGGckkkssmFFmmss..........", ".........ssmtFtsskkkkcGGGGGWWWWWGGGWWWWWWWWWWWWWWWWWWWWGGGGGGcccGGGGGcccccGGGWWWWWWWWWWWWWWWWWGGGGGGGGGGGGkkkkssmtttmss.........", "........ssmmFFssskskkcGGGGGGGGWWWWWWWWWWWWWWWWWWWWWWWWGGGcGGGGGcGGGGGGGcGGGGWWWWWWWWWWWWWWWWWWWWWWWGGGccGckkkkkmsmFtmmss........", "........smmtFtssssskkcccGGGGGGGWWWWWWWWWWWWWWWWWWWWWGGGGcGGGGGGcGGGccccGGGGGWWWWWWWWWWWWWWWWWWGGGGGGGGccGGckkkkmmsFtmmss........", ".......ssmttFmsskkkccGGGGGGcGGGGWWWWWWWWWWWWWWWWWWWGGGGGGGGGGGGcGGGGGccGGGGGWWWWWWWWWWWWWWWGGGGGGGGGGcGGGGGckkksmmttttmss.......", ".......ssttFmsskkkkcGGGGGGGGGGGGGGWWWWWWWWWWWWWWWWGGcGGGGWWGGGcccGGGGGGcccGGGGWWWWWWWWWWWWGGGGGGGGGGGGGGWGGGckkkmmsmtttmss......", ".......smttFmmskkkkcGGGWWGGGGGGGGGGWWWWWWWWWWWWWWWGGGGGGWWWGGGcccGGGGGGcccGGGGGWWWWWWWWWWGGGGGGGGGGGGWWWWWGGckkksmmmtFtmss......", "......smmtFtsmkkkccGGGGWWWWWGGGGGGGGWWWWWWWWWWWWWWWGGGWWWWWGGGcccGGWWGGGGGGGccGWWWWWWWWWWGGGGWWWGGGGGGGGWWGGGkkkksmmmtFmsss.....", "......smmFFtmmkkkGGGGGGGGWWWWWWGGGGWWWWWWWWWgWWWWWWWWWWWWGGGGGcccGGWWGGGGGGGGGGWWWWWWWWWWGGWWWGGGGGGGGGGWWWGGckkksmmmmFtmsss....", ".....smmtttmmskkcGGGGGGGGGGWWWWWWWWWWWWWWWWWcgWWWWWWWWWWGGGGGGccccGGWWGGGGGWWWWWWWWWWWWWWWWWWWGGGccGGWWWWWWGGckkkkmmmmttmsss....", ".....smmttmmmskccGGWGGGGGGGGGWWWWWWWWWWWWWWgccgWWWWWWWWWGGGGGGccccGGGWWGGGWWWWWWWWWWWWWWWWWWWWGGGGGGGWWWWWGGGcckkkmmmstttsss....", "....ssmmttmmskkccGGWWWGGGGGGGWWWWWWWWWWWWWWgcccWWWWWWWWWWGGGGGccccGGGWWWWWWWWWWWWgggWWWWWWWWWWGGGGGGWWWWWGGGGGckkksmmmmFtsss....", "....ssmtFtmmskkccGGWWWWWWWWWWWWWWWWWWWWWWWWgcccgWWWWWWWWWWWGGcccccGGGWWWWWWWWWWWgcccWWWWWWWWWWGGGWWWWWWWGGGGGGGckksmmmmFtsss....", "....ssmtFmmmskkccGGGGGWWWWWWWWWWWWWWWWWWWWWgccccgWWWWWWGGGGGGcccGGGGGWWWWWWWWWWgcccgWWWWWWWWWWWWWWWWWGGGGGGGGGGckksmmmmttsss....", "...ssmmttmmmkkcccGGGGGGGGGGWWWWWWWWWWWWWWWWgcccccWWWWWWWGGGGGcccGGGGGGGWWWWWWWgccccgWWWWWWWWWWWWWWWWGGGGGGGGWGGGckksmmmttsss....", "...ssmtttmmmkkkccGGGGGGGGGGWWWWWWWWWWWWWWWWgcccccgWWWWWWWWWccccccGGGGWWWWWWWWgcccccWWWWWWWWWWWWWWWWWWWGGGWWWWGGGckksmmmttmsss...", "..sssmtttmmmkkkcGGGGGGGGGWWWWWGWWWWWWWWWWWWgccccccgWWWWWWWWGgcccccGGWWWWWWWWgccccccWWWWWWWWWWWWWWWWWWWWWWWWWWGGGckksmmmttmsss...", "..ssmmtttmmmkkccGGGGGGGGGGGGGGGGGGGWWWWWWWWWccccccccWWWWWWWWWgggggWWWWWWWWWWccccccgWWWWWWWWWWWWWWWWWWWWWGGWWWWGGckkkmmmttmsss...", "..sssmtttmmmkkkcGGGGGGGGGGGGGGccGGGGWWWWWWWWccccccccgWWWWWWWWWWWWWWWWWWWWWgcccccccWWWWWWWWWGGGGGWWWWWWWWGGGGWWGGckkkmmmttmsss...", "..sssmtttmmmkkkcGGGGGGGGGGGGGcccGGGGWWWWWWWWgcccccccggWWWWWWWWWWWWWWWWWWWWcccccccgWWWWWWWWWGGGGGGGGGGGWWGGGGWWGGckkkmmmmtmsss...", "..sssttttmmmkkkccGGGGGcccccccGGGGGWWGGGWWWWWWgccccccccgWWWWWWWWWWWWWWWWWggcccccccWWWWWWWWWWGGGGGGGGGGGGGGGGGGGGGcckkmmmmtmsss...", "..ssmtttmmmmkkkccGGGccccGGGGGGGWWWWGGccGGWWWWWGccccccccgWWWWWWWWWWWWWWWgcccccccgWWWWWWWWGGGGWWGGGGGcccccGGGGGGGGcckkmmmmtmsss...", "..ssmtttmmmmkkkcGGGGGGGGGGGGGWWWWWGGGccGGWWWWWWcccccccccgWWWWWWWWWWWWWgccccccccWWWWWWWWWGGGGGGWGGGGGccccGGGGGGGGckkkmmmmttmss...", "...sstttmmmmkkkcGGGGGGGGGGGGWWWWGGGGccccGGWWWWWGGccccccccWWccgggggcgggcccccccggWWWWWWWGGGcccGGGGGGGGGGGGGGGGGGGGcckkmmmmtmsss...", "...ssmttmmmmkkkGGGWWWWWWWWWWWWWGGGGGGGGGGGWWWWWWGGccccccccggccccccgWgcccGGccGgWWWWWWWGccGccGGGGGGGWWWWWGWWGGGGGGcckkmmmmttmss...", "...ssmttmmmmkkkcGGGWWWWWWWWWWWGGGGGGGGGGGWWWWWWWWGGcccccccggccccccggccccGGccGWWWWWWWWGcccccGGGGGGWWWGGGGGGGGWGGGGckkmmmmmmmss...", "...ssmttmmmmkkkcGGGGGGGGGWWWWGGGGGWGGGGGWWWWWWWWWGGGGGGGGccWgcccccgcccGGGGcGWWWWWWWWWGGGGGGGGGGGWWWGGGGGGGGGGWGGGckkmmmmmtmss...", "...ssmttmmmmkkcGGGGcccGGGGWWWGGGGWGGGcGGWWWWWWWWWGGWGGGGGccgWgcccgggGccGGGGGWWWWWWWWWGGGGGGGGGGGGGWWWWWWWWWWWWGGGckkmmmmmtmss...", "....smttmmmmkkcGGGGGccGGGWWWWGGGWWGGGGGGWWWWWWWWWGGWWWGGGGcgWgccgWggGcGGGGGWWWWWWWWWWGGGGGWGGGGGGGGWWWWWWWWGGGGGGckkmmmmttms....", "....smttmmmmkkGGGWWGGGGWWWWWWWWWWGGGGGGGWWWWWWWWWWWWWWGGGGcGWWccWWcccGGGGGWWWWWWWWWWWWGGGGGGGWGGGGGGWWWWWWGGGGGGcckkmmmmttms....", "....smttmmmmkkGGGWWWWWWWWWWWWWGGGGGGGGGWWWWWWWWWWWWWWWGGGGccGWgcWgcccGGGGGWWWWWWWWWWWWWGGccGGGGGWGGGGWWWWGGGcGGGGckkmmmmtmss....", "....smttmmmmkkGGGWWWWWWWWWWWGGGGGGGccGGWWWWWWWWWWWWWWWGGGGGcGWggWgccccccGGGWWWWWWWWWWWWGGGGGGGGGWWWGGWWWGGGGGGGGGckkmmmmtmss....", "....smmtmmmmkkcGGWWWWWWWGGGGGGGGGGGGGGGWWWWWWWWWWWWWGGGGGGGccgWgWgccccccGGGWWWWWWWWWWWWGGGGGGGGGWWWWWWWWGGGGGGWGGckkmmmmtmss....", "....smmFmmmmkkcGGWWWWWGWWWWGGGWWWWGGGGGWWWWWWWWWWWWWGGGccGGccgWWgcccccccGGGGWWWWWWWWWWWWGGGGGGGWWWWWWWWWWWGGWWWGGckkmmmmtmss....", "....smmtmmmmskcGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGcccGccccWWccccccccGGGGWWWWWWWWWWWWGGGGGGGGGGGWWWWWWWWWWWGGGckkmmmmtmss....", "....smmtmmmmskcGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGGccGccccgWccccccccGGGGWWWWWWWWWWWWGGGGGGGGGGGWWWWWWWWWWGGGGckkmmmmmmss....", "....ssmtmmmmskkcGGGWWWWWGGGGWWWWWWWWWWWWWWWWWWWWWWGGGGGcccccccggccccccccGGGGWWWWWWWWWWWWWWWWWGGGGGGWWWWWWWWWGGGGcckkmmmmtsss....", "....ssmtmmmmskkscGGWWGGGcGGGWWWWWWWWWWWWWWWWWWWWWWGGGGGccccccccgccccccccGGGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGcckksmmmmtsss....", "....ssmtmmmmskkkcGGWWGcccGGGWWWWWWWWWWWWWWWWWWWWWGGGGGGccccccccgcccccccGGGGGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGcckksmmmmtsss....", "....ssmttmmmmskkcGGGGGcGGGGWWWWWWWWWWWWWWWWWWWWWWWGGGGGccccccccgcccccccGGGGGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGcckksmmmmtsss....", ".....smttmmmmskkkGGGGccGGGGGGGGWWWWWWWWWWWWWWWWWWGGGGGGGcccccccgcccccccGGGGGGWWWWWWWWWWWWWWWWWWWWGGGWWWWWWGGGGGckkkmmmmttss.....", ".....smttmmmmskkkcccccccGGGGGGGGWWWWWWWWWWWWWWWWWGGGGGGGcccccccgcccccccGGGGGGWWWWWWWWWWWWWWWWWWWGGGGGWWWWWGGGccckkkmmmmttss.....", ".....smttmmmmskkkcccGGGGGGGGGGGWWWWWWWWWWWWWWWWWWWGGGGGGcccccccgcccccccGGGGGGGWWWWWWWWWWWWWWWWWWGGccGGGGWWWGGccckksmmmmttss.....", ".....ssttmmmmmkkkccGGGGGGGGGGWWWWWWWWWWWWWWWWWWWWWGGGGGccccccccgcccccccGGGGGGGWWWWWWWWWWWWWWWWWWGGGGGcGGGGWWWGGkkksmmmmttss.....", ".....ssttmmmmmkkkccGGGGGWWGWWWWWWWWWWWWWWWWWWWWWWWGGGGGcccccccggcccccccGGGGGGGWWWWWWWWWWWWWWWWWWGGGGGGGcGGGWWGGkkksmmmmttss.....", ".....ssttmmmmmkkscccccGGGWWWWWWWWWWWWWWWWWWWWWWWWWGGGGccccccccgWggccccccGGGGGGWWWWWWWWWWWWWWWWWWWWWGGGGGGGGGGGGskkmmmmmtts......", ".....sstttmmmmkkssGGGGGGWWWWWWWWWWWWWWWWWWWWWWWWWWGGGGccccccggWWWWgggcccccGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGcccccskkmmmmmtts......", ".....sstttmmmmkkssGGGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGccccccggWWWWWWWgcccccGWWWWWWWWWWWWWWWWWWWWWWWWWGGGGGGccccckkkmmmmmttss.....", "......smttmmmmskksmGGWWWWGGGWWWWWWWWWWWWWWWWWWWWWWWWcccccgWWWWWWWWWWWWgcccGWWWWWWWWWWWWWWWWWWWWWWWWGGGGGGGccccskkkmmmmmtts......", "......sstttmmmskkssGGGGGGGGGGWWWWWWWWWWWWWWWWWWWWWWWcccggWWWWWWWWWWWWWWWWgWWWWWWWWWWWWWWWWWWWWWWWWWGGGGGGGGGGcskkkmmmmttms......", "......sstttmmmskkssccccccccGGGWWWWWWWWWWWWWWWWWWWWWWccggWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGWWGGGsskksmmmmttms......", ".......smttmmmmkkkscGGGGGGGGGWWWWWWWWWWWWWWWWWWWWWWWggWWWWWWWWWGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGsskksmmmmttms......", ".......sstttmmmskksmGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGccccGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGGccskksmmmmttmm.......", ".......sstttmmmmkkksGGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGccccGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGGccckkksmmmmttmm.......", "........smttmmmmskkscGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGGGcccGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGGccskksmmmmmttms.......", "........ssmtmmmmskksscGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGGcGGGGGWWWWWWWWWWWWWWWWWWWWWWWWWGGGGGGWWWWGGcckkksmmmmttmms.......", ".........smmtmmmmkkkksGGWWWGGGGGGGGGGGWWWWWWWWWWWWWWWWWWWWWWGGccGGGGWWWWWWWWWWWWWWWWWWWWWWWWWGGGcccGGGWWGGckkkkmmmmmttms........", ".........ssmtmmmmkkkksGGWWWGGGGGGGGGGGGWWWWWWWWWWWWWWWWWWWWGGGccGGGWWWWWWWWWWWWWWWWWWWWWWWWWWWGGGcccGGGGGGckkkkmmmmmtmms........", "..........ssttmmmskkkkGGGGGGGGGGGGGGGGGWWWWWWWWWWWWWWWWWWWGGGGGcGGGWWGGGGGWWWWWWWWWWWWWWWWWWWWGGGGGccGGGcckkkksmmmmtmms.........", "..........ssmtmmmmskkkkcGGGGGGGGGGGGGGGWWWWWWWWWWWWWWWWWWWGGGGccGGGGGGGGGGWWWWWWWWWWWWWWWWWWWWWWWGGGGcccckkkkkmmmmttmss.........", "...........smmtmmmskkkkccGGGGGGGGGGGGGWWWWWWWWWWWWWWWWWWWWGGGGccGGGGGGGGGGWWWWWWWWWWWWWGGWWWWWWWWWGGGccccsskksmmmmttmss.........", "............smttmmmkkkkkcccGGGGGGGGWWWWWWWWWWWWWWWWWWWGGGWWGGGGGGGGGccGGGGWWWWWWWWWWWWWGGGGWWWWWWGGGccccssskkmmmmttmss..........", "............ssmttmmskksskccGGGGGWWWWWWWWWWWWWWWWWWWWWWGGGGGWWGGGcccGGGGGGGWWWWWWWWWWWWWWGGGGGGGWWGGGcccksskksmmmttmss...........", ".............smttmmmsksskkcGGGGGGWWWWWWWWWWWWWWWWWWWWWGGGGGGWGGGcccGGGGGGGWWWWWWWWWWWWWWWGGGGGGWWGGGcckkssksmmmtFtss............", "..............smttmmmskskkcGGGGGGWWWWWWWWWGGGWWWWWWWWWGGGGGGGGGGcccccGGGGGWWWWWWWWWWWWWWWWGGGGWWWGGGckkskksmmmmtmsss............", "..............ssmttmmmkkskkGGGWWWWWWWWWWGGGGGWWWWWWWWWGGGGGGGGGcccGGGGGWWWWWWWWWWWGGGGWWWWWGGGWWGGGcckksksmmmtFmss..............", "..............sssttmmmskskkcGGWWWWWWWWGGGGGGWWWWWWWWWWWGGGGGGGcccGGGGGWWWWWWWWWWWWGGGGGWWWWGGWWWGGGcckssksmmmFFmss..............", "................ssmtmmmskkkkcGGWWWGGGGGGGcGGWWWWWWWWWWWWGGGGGGGcGGGGWWWWWWWWWWWWWGGGGcGGWWWWWWWGGGcckssksmmmmtmss...............", ".................sstttsskkkkkGGGGGGGGGGcccGGWWWWWWWWWWWWWGGGGGGGGGGGGGGGWWWWWWWWGGGGGGGGGWWWWWGGGcckksksmmmFFmss................", ".................ssmtFsskkkkkGGGGGGGGGGcccGGWWWWWWWWWWWWWGGGGGGGcGGGGGGGWWWWWWWWGGGGGGcGGWWWWWGGccksskksmmtFtsss................", "..................ssmttmskkkkccGGGGGGGGGccGGGGWWWWWWWWWWWWWWGGGGccGGGGGWWWWWWWWGGGGGGGcGGWWGGGGGckkssksmmtFmsss.................", "...................sssFFsskkkkkccGGGGGGGccGGGGGWWWWWWWWWWWWWGGGGcGGGGWWWWWWWWGGGGGGGGcGGGWGGGGckkksskssstFsss...................", "...................ssstFmsskkskkcGGGGGGGcccGGGGWWWWWWWWWWWWWGGGGcGGGWWWWWWWWGGGGGGGGccGGGGGGGckkssskkssmFtsms...................", ".....................ssmFtmkkkkskkcGGGGGcccccGGGWWWWWWWWWWWWGGGGcGGGGGWWWWWWGGcccGGGccccGGGcckkssskksmttmsss....................", "......................sstFFmskkksskccGGcccccGGGWWWWWWWWWWGGGGGGGcGGGGGGGGWWWGGcccGGGGcccccckkksskkksmFFtss......................", ".......................ssmtFtsskksskcccccccccGGWWWWGGGGWWWGGGGGccGGGGGGGGGWWGGcccGGGGGccckskksskksstFtmsss......................", "........................ssmFFmsskssskccccccccGGWWWWGGGGGWWGGGGcccGGGGGGGGGWWGGccccGGcccckkkksskkssmFFmsss.......................", "........................sssmtFttsksmmskccccccGGGGGGGcccGGWWGGGccccGGGGccGGWWGGcccccccckkkksskkksttFtmsss........................", "..........................sssmtttmkkssssskcccccGGcccccccGGWWGGGcccGGGGccGGGWGGGGcccckkskssskkstttttsss..........................", "...........................sssmtttsskssssskccccccccccccccGGGGGGcccGGGGccGGGGGGGGcccckssssskksmFFtmmss...........................", ".............................ssmmtFtmkkssskkcccccmcccccccGGGGGcsccccccccccGGGGccckkssssssksttFtmmsss............................", "..............................sssmttFtmkkkkkkkkkkkccccccccccccssccccccccccccckkkssssmskksmttttmmss..............................", "................................ssmmtFtmsskkkkkssskkkkccccccccsmscccccccckkkkkksssssssssmtttmmmss...............................", "..................................smmmtFFFmskkkssssssskkkccccsssssscckkkssssssssskkksmtFFtmssss.................................", "....................................ssmmttFFttmsskkskssssskkkkkkkkkssssssssskkkkssmtFFttmmssss..................................", ".....................................ssmmmttFFFtmsskkkkssssssssssssssssssskkkkssmttFFtmmmsss....................................", ".......................................sssmmmmmtFFFtmmsskkkkkkkkkkkkkkkkksssmttFFtttmmssss......................................", "..........................................sssmmmmttttFFFttmmsssssssssmttFFtttttmmmmmsss.........................................", "............................................sssmmmmmtttFFFFFFtttmmmmtFFFFFttmmmmmmssss..........................................", "...............................................sssssmmmmmmmmmttttmmmmmmmmmmsmsssss..............................................", "....................................................ssssssssmsmmmmmmmmsssssss...................................................", ".........................................................ssssssssssssssss.......................................................", "................................................................................................................................", "................................................................................................................................", }; /* =============================== cut here ============================= */ /* ============================== end brainmri128.h ========================= */ /* =============================== cut here ============================= */