~ubuntu-branches/ubuntu/utopic/mricron/utopic

function GetValX (var lnSubj, lnFactors: integer; var lSymptomRA: singleP; var lImageNames: TStrings; var lCrit: integer; {lBinomial : boolean;} var lPredictorList: TStringList):boolean;

function ThreshMap(lThresh: single; lVolVox: integer;lOutImg: singleP): integer;

function FirthNPMAnalyze (var lImages: TStrings; var lPredictorList: TStringList; var lMaskHdr: TMRIcroHdr; lnCond,lnCrit: integer; var lSymptomRA: SingleP; var lOutName: string): boolean;

100

procedure InitPermute (lnSubj, lnPermute: integer; var lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP; var lRanOrderp: pointer; var lRanOrder: Doublep0);

101

function reportPermute (lLabel:string; lnPermute: integer; var lPermuteMaxZ, lPermuteMinZ: singleP): double;

102

procedure FreePermute (lnPermute: integer; var lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP;var lRanOrderp: pointer);

103

procedure FormClose(Sender: TObject; var CloseAction: TCloseAction);

104

function SaveHdrName (lCaption: string; var lFilename: string): boolean;

105

procedure StrToMemo (lStr: string);

106

procedure NPMclick(Sender: TObject);

107

function OpenDialogExecute (lCaption: string;lAllowMultiSelect,lForceMultiSelect: boolean; lFilter: string): boolean;//; lAllowMultiSelect: boolean): boolean;

108

function NPMAnalyze (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lMaskVoxels,lnGroup1: integer): boolean;

109

function NPMAnalyzePaired (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lMaskVoxels: integer): boolean;

110

function NPMzscore (var lImages: TStrings; var lMnHdr,lStDevHdr: TMRIcroHdr): boolean;

111

procedure FormCreate(Sender: TObject);

112

function ReportDescriptives (var RA: SingleP; n: integer): boolean;

113

function MakeSubtract (lPosName,lNegName: string): boolean;

114

function MakeMean (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lBinarize,lVariance: boolean): boolean;

115

function Balance (var lImageName,lMaskName: String; lMethod: integer{lInflection: boolean}): boolean;

116

procedure LesionBtnClick(Sender: TObject);

117

procedure Copy1Click(Sender: TObject);

118

procedure testmenuclick(Sender: TObject);

119

procedure radiomenuclick(Sender: TObject);

120

procedure ReadIniFile;

121

procedure WriteIniFile;

122

function reportBonferroni(lLabel: string; lnTests: integer): double;

123

function reportFDR (lLabel:string; lnVox, lnTests: integer; var lData: SingleP): double;

124

procedure Makemeanimage1Click(Sender: TObject);

125

procedure Exit1Click(Sender: TObject);

126

procedure Balance1Click(Sender: TObject);

127

procedure niiniigz1Click(Sender: TObject);

128

procedure Variance1Click(Sender: TObject);

129

procedure ZtoP1Click(Sender: TObject);

130

procedure About1Click(Sender: TObject);

131

procedure Design1Click(Sender: TObject);

132

procedure PhysiologicalArtifactCorrection1Click(Sender: TObject);

133

procedure DualImageCorrelation1Click(Sender: TObject);

134

procedure FormShow(Sender: TObject);

135

procedure PairedTMenuClick(Sender: TObject);

136

procedure SingleSubjectZScores1Click(Sender: TObject);

137

procedure MultipleRegressClick(Sender: TObject);

138

function ReadPermute: integer;

139

procedure NPMmsg( lStr: string);

140

procedure NPMmsgClear;

141

procedure MsgSave(lFilename: string);

142

procedure SingleRegressClick(Sender: TObject);

143

procedure AssociatevalfileswithNPM1Click(Sender: TObject);

144

procedure threadChange(Sender: TObject);

145

procedure Countlesionoverlaps1Click(Sender: TObject);

146

procedure PenalizedLogisticRegerssion1Click(Sender: TObject);

147

procedure ComputeIntersectionandUnion1Click(Sender: TObject);

148

procedure ROCbinomialdeficit1Click(Sender: TObject);

149

procedure ROCcontinuousdeficit1Click(Sender: TObject);

150

procedure ThreadDone(Sender: TObject);

151

procedure NPMmsgAppend( lStr: string);

152

procedure ROIanalysis1Click(Sender: TObject);

153

procedure Masked1Click(Sender: TObject);

154

procedure Binarizeimages1Click(Sender: TObject);

155

procedure Resliceimagetoneworientationandboundingbox1Click(

156

Sender: TObject);

157

procedure Setnonseroto1001Click(Sender: TObject);

158

procedure Savetext1Click(Sender: TObject);

159

procedure AnaCOMmenuClick(Sender: TObject);

160

procedure MonteCarloSimulation1Click(Sender: TObject);

161

procedure Subtract1Click(Sender: TObject);

162

procedure LogPtoZ1Click(Sender: TObject);

163

private

164

{ Private declarations }

165

public

166

{ Public declarations }

167

end;

168

169

var

170

MainForm: TMainForm;

171

implementation

172

173

uses filename,prefs,roc,hdr,regression,valformat {$IFDEF SPREADSHEET} ,design,spread{$ENDIF}

174

{$IFNDEF UNIX},ActiveX {$ENDIF};

175

{$IFNDEF FPC}

176

{$R *.DFM}

177

178

{$ENDIF}

179

const

180

kVers : string = 'Chris Rorden''s NPM :: '+kMRIcronVers;

181

var

182

gNULP: boolean = true;

183

gROI : boolean = false;

184

function TMainForm.GetKVers: string;

185

begin

186

result := kVers +'; Threads used = '+inttostr(gnCPUThreads );

187

end;

188

189

190

procedure TMainForm.NPMmsgAppend( lStr: string);

191

var

192

lOutname: string;

193

f: TextFile;

194

begin

195

MainForm.Memo1.Lines.add(lStr);

196

lOutname:='c:\dx.txt';

197

if fileexists(lOutname) then

198

begin { open a text file }

199

AssignFile(f, lOutname);

200

Append(f);

201

Writeln(f, lStr);

202

Flush(f); { ensures that the text was actually written to file }

203

{ insert code here that would require a Flush before closing the file }

204

CloseFile(f);

205

end;

206

end;

207

208

procedure TMainForm.NPMmsg( lStr: string);

209

begin

210

MainForm.Memo1.Lines.add(lStr);

211

end;

212

213

procedure Msg(lStr: string);

214

begin

215

MainForm.NPMmsg(lStr);

216

end;

217

218

procedure MsgClear;

219

begin

220

MainForm.Memo1.Lines.Clear;

221

end;

222

223

procedure TMainForm.NPMmsgClear;

224

begin

225

MsgClear;

226

end;

227

228

229

procedure TMainForm.MsgSave(lFilename: string);

230

begin

231

MainForm.Memo1.Lines.SaveToFile(lFilename);

232

end;

233

234

procedure TMainForm.ThreadDone(Sender: TObject);

235

begin

236

Dec(gThreadsRunning);

237

end;

238

239

procedure InitRA (lnPermute: integer; var lRA: singleP);

240

var

241

lInc: integer;

242

begin

243

getmem(lRA,lnPermute* sizeof(single));

244

for lInc := 1 to lnPermute do

245

lRA^[lInc] := 0;

246

end;

247

248

procedure TMainForm.InitPermute (lnSubj, lnPermute: integer; var lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP; var lRanOrderp: pointer; var lRanOrder: Doublep0);

249

begin

250

if (lnPermute < 2) then

251

exit;

252

InitRA(lnPermute,lPermuteMaxT);

253

InitRA(lnPermute,lPermuteMinT);

254

InitRA(lnPermute,lPermuteMaxBM);

255

InitRA(lnPermute,lPermuteMinBM);

256

createArray64(lRanOrderp,lRanOrder,lnSubj);

257

end; //init permute

258

259

procedure sort (lo, up: integer; var r:SingleP);

260

//62ms Shell Sort http://www.dcc.uchile.cl/~rbaeza/handbook/algs/4/414.sort.p.html

261

label 999;

262

var d, i, j : integer;

263

tempr : single;

264

begin

265

d := up-lo+1;

266

while d>1 do begin

267

if d<5 then

268

d := 1

269

else

270

d := trunc( 0.45454*d );

271

//Do linear insertion sort in steps size d

272

for i:=up-d downto lo do begin

273

tempr := r^[i];

274

j := i+d;

275

while j <= up do

276

if tempr > r^[j] then begin

277

r^[j-d] := r^[j];

278

j := j+d

279

end

280

else goto 999; //break

281

999:

282

r^[j-d] := tempr

283

end //for

284

end //while

285

end; //proc Sort

286

287

function IndexPct(lnPermute: integer; lPct: single; lTop: boolean): integer;

288

begin

289

result := round(lnPermute * lPct);

290

if lTop then

291

result := (lnPermute - result)+1;

292

if (result < 1) then

293

result := 1;

294

if (result > lnPermute) then

295

result := lnPermute;

296

end;

297

298

function TMainForm.reportBonferroni(lLabel: string; lnTests: integer): double; //returns 5% Z score

299

begin

300

if lnTests < 1 then exit;

301

result := pNormalInv(0.05/lnTests);

302

msg(inttostr(lnTests)+' test '+lLabel+' Bonferroni FWE Z '+

303

'0.050='+realtostr(result,3)+

304

', 0.025='+realtostr(pNormalInv(0.025/lnTests),3)+

305

', 0.01='+realtostr(pNormalInv(0.01/lnTests),3));

306

end;

307

308

function TMainForm.reportFDR (lLabel:string; lnVox, lnTests: integer; var lData: SingleP): double;

309

var

310

lC,lN: integer;

311

lPs: SingleP;

312

lFDR05r, lFDR01r,lFDR05p, lFDR01p,lMin,lMax : double;

313

begin

314

result := 10000;

315

if (lnTests < 1) or (lnVox < 1) then

316

exit;

317

GetMem(lPs,lnTests*sizeof(single));

318

for lC := 1 to lnTests do

319

lPs^[lC] := 0;

320

lN := 0;

321

lMin := 0;

322

lMax := 0;

323

for lC := 1 to lnVox do begin

324

if lData^[lC] <> 0 then begin

325

inc(lN);

326

if lData^[lC] > lMax then lMax := lData^[lC]

327

else if lData^[lC] < lMin then lMin := lData^[lC];

328

if lN <= lnTests then

329

lPs^[lN] := pNormal(lData^[lC]);

330

end;

331

end;

332

EstimateFDR2(lnTests, lPs, lFDR05p, lFDR01p,lFDR05r, lFDR01r);

333

msg(lLabel+' Range '

334

+realtostr(lMin,3)+

335

'...'+realtostr(lMax,3));

336

{Msg(lLabel+' Range '

337

+realtostr(pNormalInv(lPs[lnTests]),3)+

338

'...'+realtostr(pNormalInv(lPs[1]),3)+

339

' '); } //we could use this and save time computing lmin/lmax, but loss in precision

340

msg(lLabel+' +FDR Z '+

341

'0.050='+realtostr(pNormalInv(lFDR05p),8)+

342

', 0.01='+realtostr(pNormalInv(lFDR01p),8)+

343

' ');

344

msg(lLabel+' -FDR Z '+

345

'0.050='+realtostr(pNormalInv(1-lFDR05r),8)+

346

', 0.01='+realtostr(pNormalInv(1-lFDR01r),8)+

347

' ');

348

result := pNormalInv(lFDR01p);

349

end;

350

351

function ReportThresh (lLabel: string; lnPermute: integer; var lRankedData: singleP;lTop:boolean): double;

352

begin

353

result := lRankedData^[IndexPct(lnPermute,0.050,lTop)];

354

msg(lLabel+': permutationFWE '+

355

//'0.500='+realtostr(lRankedData[IndexPct(lnPermute,0.500,lTop)],3)+

356

', 0.050='+realtostr({lRankedData^[IndexPct(lnPermute,0.050,lTop)]} result,8)+

357

', 0.025='+realtostr(lRankedData^[IndexPct(lnPermute,0.025,lTop)],8)+

358

', 0.01='+realtostr(lRankedData^[IndexPct(lnPermute,0.010,lTop)],8)+

359

' ');

360

end;

361

362

function TMainForm.reportPermute (lLabel:string; lnPermute: integer; var lPermuteMaxZ, lPermuteMinZ: singleP): double;

363

begin

364

result := 0;

365

if (lnPermute < 2) then

366

exit;

367

sort (1, lnPermute,lPermuteMaxZ);

368

result := ReportThresh(lLabel+'+',lnPermute,lPermuteMaxZ,true);

369

sort (1, lnPermute,lPermuteMinZ);

370

ReportThresh(lLabel+'-',lnPermute,lPermuteMinZ,false);

371

//for lPos := 1 to lnPermute do

372

// msg(inttostr(lPos)+', '+realtostr(lPermuteMinZ[lPos],4));

373

374

end;

375

376

procedure TMainForm.FreePermute (lnPermute: integer; var lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP;var lRanOrderp: pointer);

377

begin

378

if (lnPermute < 2) then

379

exit;

380

Freemem(lRanOrderp);

381

Freemem(lPermuteMaxT);

382

Freemem(lPermuteMinT);

383

Freemem(lPermuteMaxBM);

384

Freemem(lPermuteMinBM);

385

end;

386

387

function TMainForm.SaveHdrName (lCaption: string; var lFilename: string): boolean;

388

begin

389

result := false;

390

SaveHdrDlg.InitialDir := lFilename;

391

SaveHdrDlg.Title := lCaption;

392

SaveHdrDlg.Filter := kAnaHdrFilter;

393

if not SaveHdrDlg.Execute then exit;

394

lFilename := SaveHdrDlg.Filename;

395

result := true;

396

end;

397

398

procedure TMainForm.FormClose(Sender: TObject; var CloseAction: TCloseAction);

399

begin

400

//showmessage('zz');

401

WriteIniFile;

402

end;

403

404

procedure MakeStatHdr (var lBGHdr,lStatHdr: TniftiHdr; lMinIntensity,lMaxIntensity,lIntent_p1,lIntent_p2,lIntent_p3: single; lIntent_code: smallint;lIntentName: string);

405

var lIntentNameLen,lPos: integer;

406

lStr: string;

407

begin

408

move(lBGHdr,lStatHdr,sizeof(TniftiHdr));

409

with lStatHdr do begin

410

magic :=kNIFTI_MAGIC_SEPARATE_HDR;

411

bitpix := 32; //32-bit real data

412

datatype := kDT_FLOAT;

413

scl_slope:= 1;

414

scl_inter:= 0;

415

glmin := round(lMinIntensity);

416

glmax := round(lMaxIntensity);

417

intent_code := lIntent_Code;// kNIFTI_INTENT_ESTIMATE;

418

intent_p1 := lIntent_p1;

419

intent_p2 := lIntent_p2;

420

intent_p3 := lIntent_p3;

421

lIntentNameLen := length(lIntentName);

422

descrip[1] := 'N';

423

descrip[2] := 'P';

424

descrip[3] := 'M';

425

if lIntent_code=kNIFTI_INTENT_TTEST then begin

426

descrip[4] := 't' ;

427

lStr := inttostr(trunc(lIntent_p1));

428

for lPos := 1 to length (lStr) do

429

descrip[4+lPos] := lStr[lPos] ;

430

end else

431

descrip[4] := 'z';

432

if lIntentNameLen > sizeof(intent_name) then

433

lIntentNameLen := sizeof(intent_name);

434

if lIntentNameLen > 0 then

435

for lPos := 1 to lIntentNameLen do

436

intent_name[lPos] := lIntentName[lPos];

437

end;

438

end;

439

440

procedure WriteThread( lnThread: integer);

441

begin

442

case lnThread of

443

2: MainForm.T2.checked := true;

444

3: MainForm.T3.checked := true;

445

4: MainForm.T4.checked := true;

446

7: MainForm.T7.checked := true;

447

8: MainForm.T8.checked := true;

448

15: MainForm.T15.checked := true;

449

16: MainForm.T16.checked := true;

450

else MainForm.T1.checked := true;

451

end;

452

gnCPUThreads := lnThread;

453

end;

454

455

function ReadThread: integer;

456

begin

457

if MainForm.T16.checked then result := 16

458

else if MainForm.T15.checked then result := 15

459

else if MainForm.T8.checked then result := 8

460

else if MainForm.T7.checked then result := 7

461

else if MainForm.T4.checked then result := 4

462

else if MainForm.T3.checked then result := 3

463

else if MainForm.T2.checked then result := 2

464

else result := 1;

465

gnCPUThreads := result;

466

end;

467

468

469

procedure WritePermute( lnPermute: integer);

470

begin

471

case lnPermute of

472

4000: MainForm.N4000.checked := true;

473

3000: MainForm.N3000.checked := true;

474

2000: MainForm.N2000.checked := true;

475

1000: MainForm.N1000.checked := true;

476

else MainForm.N0.checked := true;

477

end;

478

end;

479

480

function TMainForm.ReadPermute: integer;

481

begin

482

if MainForm.N4000.checked then result := 4000

483

else if MainForm.N3000.checked then result := 3000

484

else if MainForm.N2000.checked then result := 2000

485

else if MainForm.N1000.checked then result := 1000

486

else result := 0;

487

end;

488

489

function TMainForm.NPMAnalyze (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lMaskVoxels,lnGroup1: integer): boolean;

490

label

491

667;

492

var

493

lOutName,lOutNameMod: string;

494

lMaskImg,lPlankImg,lOutImgMn,lOutImgBM,lOutImgT,lDummy: SingleP;

495

lTotalMemory: double; //not integer - limit for 32bit int is 2Gb

496

lPlank,lVolVox,lPos,lMinMask,lMaxMask,lnPlanks,lVoxPerPlank,

497

lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lThreadStart,lThreadEnd,lThreadInc: integer;

498

lT, lSum, lMn: double;

499

lStatHdr: TNIfTIhdr;

500

lFdata: file;

501

lThread,lnPermute: integer;

502

lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP;

503

lRanOrderp: pointer;

504

lRanOrder: Doublep0;

505

lttest,lBM: boolean;

506

begin

507

lttest:= ttestmenu.checked;

508

lBM := BMmenu.checked;

509

510

lnPermute := ReadPermute;

511

//lnPermute := 100;

512

msg('Permutations = ' +IntToStr(lnPermute));

513

lOutName := lMaskHdr.ImgFileName;

514

if not SaveHdrName ('Statistical Map', lOutName) then exit;

515

msg('Analysis began = ' +TimeToStr(Now));

516

lTotalMemory := 0;

517

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

518

if (lVolVox < 1) then goto 667;

519

//load mask

520

getmem(lMaskImg,lVolVox*sizeof(single));

521

if not LoadImg(lMaskHdr.ImgFileName, lMaskImg, 1, lVolVox,round(gOffsetRA[0]),1,lMaskHdr.NIFTIhdr.datatype,lVolVox) then begin

522

msg('Unable to load mask ' +lMaskHdr.ImgFileName);

523

goto 667;

524

end;

525

//next find start and end of mask

526

lPos := 0;

527

repeat

528

inc(lPos);

529

until (lMaskImg^[lPos] > 0) or (lPos = lVolVox);

530

lMinMask := lPos;

531

lPos := lVolVox+1;

532

repeat

533

dec(lPos);

534

until (lMaskImg^[lPos] > 0) or (lPos = 1);

535

lMaxMask := lPos;

536

if lMaxMask = 1 then begin

537

msg('Mask appears empty' +lMaskHdr.ImgFileName);

538

goto 667;

539

end;

540

msg('Mask has voxels from '+inttostr(lMinMask)+'..'+inttostr(lMaxMask));

541

lVoxPerPlank := kPlankSz div lImages.Count div sizeof(single) ;

542

if (lVoxPerPlank = 0) then goto 667; //no data

543

lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;

544

if (lTotalMemory = 0) then goto 667; //no data

545

lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;

546

msg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));

547

msg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));

548

getmem(lPlankImg,kPlankSz);

549

lStartVox := lMinMask;

550

lEndVox := lMinMask-1;

551

for lPos := 1 to lImages.Count do

552

if gScaleRA[lPos] = 0 then

553

gScaleRA[lPos] := 1;

554

getmem(lOutImgMn,lVolVox* sizeof(single));

555

getmem(lOutImgBM,lVolVox* sizeof(single));

556

getmem(lOutImgT,lVolVox* sizeof(single));

557

InitPermute (lImages.Count, lnPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);

558

for lPos := 1 to lVolVox do begin

559

lOutImgMn^[lPos] := 0;

560

lOutImgBM^[lPos] := 0;

561

lOutImgT^[lPos] := 0;

562

end;

563

ClearThreadData(gnCPUThreads,lnPermute);

564

for lPlank := 1 to lnPlanks do begin

565

msg('Computing plank = ' +Inttostr(lPlank));

566

Refresh;

567

Application.processmessages;

568

lEndVox := lEndVox + lVoxPerPlank;

569

if lEndVox > lMaxMask then begin

570

lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);

571

lEndVox := lMaxMask;

572

end;

573

lPlankImgPos := 1;

574

for lPos := 1 to lImages.Count do begin

575

if not LoadImg(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then

576

goto 667;

577

lPlankImgPos := lPlankImgPos + lVoxPerPlank;

578

end;//for each image

579

//threading start

580

lThreadStart := 1;

581

lThreadInc := lVoxPerPlank div gnCPUThreads;

582

lThreadEnd := lThreadInc;

583

Application.processmessages;

584

for lThread := 1 to gnCPUThreads do begin

585

if lThread = gnCPUThreads then

586

lThreadEnd := lVoxPerPlank; //avoid integer rounding error

587

with TNNStat.Create (ProgressBar1,lttest,lBM,0, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,lnGroup1, lMaskImg,lPlankImg,lOutImgMn,lOutImgBM,lOutImgT,lDummy) do

588

{$IFDEF FPC} OnTerminate := @ThreadDone; {$ELSE}OnTerminate := ThreadDone;{$ENDIF}

589

inc(gThreadsRunning);

590

lThreadStart := lThreadEnd + 1;

591

lThreadEnd :=lThreadEnd + lThreadInc;

592

end; //for each thread

593

repeat

594

Application.processmessages;

595

until gThreadsRunning = 0;

596

Application.processmessages;

597

//threading end

598

lStartVox := lEndVox + 1;

599

end;

600

lnVoxTested := SumThreadData(gnCPUThreads,lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);

601

//next report findings

602

msg('Voxels tested = ' +Inttostr(lnVoxTested));

603

reportBonferroni('Std',lnVoxTested);

604

//next: save data

605

(*savedata *)

606

MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);

607

//save mean

608

lOutNameMod := ChangeFilePostfixExt(lOutName,'Mean','.hdr');

609

if lnVoxTested > 1 then

610

611

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgMn,1);

612

MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );

613

614

if (lttest) and (lnVoxTestED > 1 ) then begin //save Ttest

615

//reportRange ('ttest', lVolVox, lnVoxTested, lOutImgT);

616

//next: convert t-scores to z scores

617

for lPos := 1 to lVolVox do

618

lOutImgT^[lPos] := TtoZ (lOutImgT^[lPos],lImages.Count-2);

619

for lPos := 1 to lnPermute do begin

620

lPermuteMaxT^[lPos] := TtoZ (lPermuteMaxT^[lPos],lImages.Count-2);

621

lPermuteMinT^[lPos] := TtoZ (lPermuteMinT^[lPos],lImages.Count-2);

622

end;

623

624

reportFDR ('ttest', lVolVox, lnVoxTested, lOutImgT);

625

reportPermute('ttest',lnPermute,lPermuteMaxT, lPermuteMinT);

626

lOutNameMod := ChangeFilePostfixExt(lOutName,'ttest','.hdr');

627

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgT,1);

628

end;

629

if (lBM) and (lnVoxTested > 1 ) then begin //save Brunner Munzel

630

reportFDR ('BM', lVolVox, lnVoxTested, lOutImgBM);

631

reportPermute('BM',lnPermute,lPermuteMaxBM, lPermuteMinBM);

632

lOutNameMod := ChangeFilePostfixExt(lOutName,'BM','.hdr');

633

634

(*reportFDR ('absT', lVolVox, lnVoxTested, lOutImgBM);

635

reportPermute('absT',lnPermute,lPermuteMaxBM, lPermuteMinBM);

636

lOutNameMod := ChangeFilePostfixExt(lOutName,'absT','.hdr');

637

638

//NIFTIhdr_SaveHdr(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr));

639

lOutNameMod := changefileext(lOutNameMod,'.img');

640

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgBM,1);

641

end;(**)

642

//next: close images

643

FreePermute (lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp);

644

freemem(lOutImgT);

645

freemem(lOutImgBM);

646

freemem(lOutImgMn);

647

//freemem(lObsp);

648

freemem(lMaskImg);

649

freemem(lPlankImg);

650

msg('Analysis finished = ' +TimeToStr(Now));

651

lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes','.txt');

652

MsgSave(lOutNameMod);

653

ProgressBar1.Position := 0;

654

exit;

655

667: //you only get here if you aborted ... free memory and report error

656

if lVolVox > 1 then freemem(lMaskImg);

657

if lTotalMemory > 1 then freemem(lPlankImg);

658

msg('Unable to complete analysis.');

659

ProgressBar1.Position := 0;

660

end;

661

662

function TMainForm.NPMAnalyzePaired (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lMaskVoxels: integer): boolean;

663

label

664

667;

665

var

666

lOutName,lOutNameMod: string;

667

lMaskImg,lPlankImg,lOutImgMn,lOutImgT,lDummy,lDummy2: SingleP;

668

lTotalMemory: double; //not integer - limit for 32bit int is 2Gb

669

lPlank,lVolVox,lPos,lMinMask,lMaxMask,lnPlanks,lVoxPerPlank,

670

lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lThreadStart,lThreadEnd,lThreadInc: integer;

671

lT, lSum, lMn: double;

672

lStatHdr: TNIfTIhdr;

673

lFdata: file;

674

lThread,lnPermute: integer;

675

lPermuteMaxT, lPermuteMinT: singleP;

676

lRanOrderp: pointer;

677

lRanOrder: Doublep0;

678

begin

679

//lnPermute := ReadPermute;

680

lnPermute := 0;//not yet

681

msg('Permutations = ' +IntToStr(lnPermute));

682

lOutName := lMaskHdr.ImgFileName;

683

if not SaveHdrName ('Statistical Map', lOutName) then exit;

684

msg('Analysis began = ' +TimeToStr(Now));

685

lTotalMemory := 0;

686

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

687

if (lVolVox < 1) then goto 667;

688

//load mask

689

getmem(lMaskImg,lVolVox*sizeof(single));

690

if not LoadImg(lMaskHdr.ImgFileName, lMaskImg, 1, lVolVox,round(gOffsetRA[0]),1,lMaskHdr.NIFTIhdr.datatype,lVolVox) then begin

691

msg('Unable to load mask ' +lMaskHdr.ImgFileName);

692

goto 667;

693

end;

694

//next find start and end of mask

695

lPos := 0;

696

repeat

697

inc(lPos);

698

until (lMaskImg^[lPos] > 0) or (lPos = lVolVox);

699

lMinMask := lPos;

700

lPos := lVolVox+1;

701

repeat

702

dec(lPos);

703

until (lMaskImg^[lPos] > 0) or (lPos = 1);

704

lMaxMask := lPos;

705

if lMaxMask = 1 then begin

706

Msg('Mask appears empty' +lMaskHdr.ImgFileName);

707

goto 667;

708

end;

709

Msg('Mask has voxels from '+inttostr(lMinMask)+'..'+inttostr(lMaxMask));

710

lVoxPerPlank := kPlankSz div lImages.Count div sizeof(single) ;

711

if (lVoxPerPlank = 0) then goto 667; //no data

712

lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;

713

if (lTotalMemory = 0) then goto 667; //no data

714

lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;

715

Msg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));

716

Msg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));

717

getmem(lPlankImg,kPlankSz);

718

lStartVox := lMinMask;

719

lEndVox := lMinMask-1;

720

for lPos := 1 to lImages.Count do

721

if gScaleRA[lPos] = 0 then

722

gScaleRA[lPos] := 1;

723

getmem(lOutImgMn,lVolVox* sizeof(single));

724

getmem(lOutImgT,lVolVox* sizeof(single));

725

//not yet InitPermute (lImages.Count, lnPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);

726

for lPos := 1 to lVolVox do begin

727

lOutImgMn^[lPos] := 0;

728

lOutImgT^[lPos] := 0;

729

end;

730

ClearThreadData(gnCPUThreads,lnPermute);

731

for lPlank := 1 to lnPlanks do begin

732

Msg('Computing plank = ' +Inttostr(lPlank));

733

Refresh;

734

Application.processmessages;

735

lEndVox := lEndVox + lVoxPerPlank;

736

if lEndVox > lMaxMask then begin

737

lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);

738

lEndVox := lMaxMask;

739

end;

740

lPlankImgPos := 1;

741

for lPos := 1 to lImages.Count do begin

742

if not LoadImg(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then

743

goto 667;

744

lPlankImgPos := lPlankImgPos + lVoxPerPlank;

745

end;//for each image

746

//threading start

747

lThreadStart := 1;

748

lThreadInc := lVoxPerPlank div gnCPUThreads;

749

lThreadEnd := lThreadInc;

750

Application.processmessages;

751

for lThread := 1 to gnCPUThreads do begin

752

if lThread = gnCPUThreads then

753

lThreadEnd := lVoxPerPlank; //avoid integer rounding error

754

with TPairedTStat.Create (ProgressBar1,false,false,0, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,666, lMaskImg,lPlankImg,lOutImgMn,lDummy2,lOutImgT,lDummy) do

755

{$IFDEF FPC} OnTerminate := @ThreadDone; {$ELSE}OnTerminate := ThreadDone;{$ENDIF}

756

inc(gThreadsRunning);

757

lThreadStart := lThreadEnd + 1;

758

lThreadEnd :=lThreadEnd + lThreadInc;

759

end; //for each thread

760

repeat

761

Application.processmessages;

762

until gThreadsRunning = 0;

763

Application.processmessages;

764

//threading end

765

lStartVox := lEndVox + 1;

766

end;

767

lnVoxTested := SumThreadDataLite(gnCPUThreads);//not yet SumThreadData(gnCPUThreads,lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);

768

//next report findings

769

Msg('Voxels tested = ' +Inttostr(lnVoxTested));

770

reportBonferroni('Std',lnVoxTested);

771

//next: save data

772

(*savedata *)

773

MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);

774

//save mean

775

lOutNameMod := ChangeFilePostfixExt(lOutName,'Mean','.hdr');

776

if lnVoxTested > 1 then

777

778

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgMn,1);

779

MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );

780

781

if (lnVoxTestED > 1 ) then begin //save Ttest

782

//next: convert t-scores to z scores

783

for lPos := 1 to lVolVox do

784

lOutImgT^[lPos] := TtoZ (lOutImgT^[lPos],(lImages.Count div 2)-1);

785

for lPos := 1 to lnPermute do begin

786

lPermuteMaxT^[lPos] := TtoZ (lPermuteMaxT^[lPos],lImages.Count-2);

787

lPermuteMinT^[lPos] := TtoZ (lPermuteMinT^[lPos],lImages.Count-2);

788

end;

789

790

reportFDR ('ttest', lVolVox, lnVoxTested, lOutImgT);

791

reportPermute('ttest',lnPermute,lPermuteMaxT, lPermuteMinT);

792

lOutNameMod := ChangeFilePostfixExt(lOutName,'ttest','.hdr');

793

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgT,1);

794

end;

795

//next: close images

796

//not yet FreePermute (lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp);

797

freemem(lOutImgT);

798

freemem(lOutImgMn);

799

//freemem(lObsp);

800

freemem(lMaskImg);

801

freemem(lPlankImg);

802

Msg('Analysis finished = ' +TimeToStr(Now));

803

lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes','.txt');

804

MsgSave(lOutNameMod);

805

ProgressBar1.Position := 0;

806

exit;

807

667: //you only get here if you aborted ... free memory and report error

808

if lVolVox > 1 then freemem(lMaskImg);

809

if lTotalMemory > 1 then freemem(lPlankImg);

810

Msg('Unable to complete analysis.');

811

ProgressBar1.Position := 0;

812

end;

813

814

815

816

function TMainForm.OpenDialogExecute (lCaption: string;lAllowMultiSelect,lForceMultiSelect: boolean; lFilter: string): boolean;//; lAllowMultiSelect: boolean): boolean;

817

var

818

lNumberofFiles: integer;

819

begin

820

OpenHdrDlg.Filter := lFilter;//kAnaHdrFilter;//lFilter;

821

OpenHdrDlg.FilterIndex := 1;

822

OpenHdrDlg.Title := lCaption;

823

if lAllowMultiSelect then

824

OpenHdrDlg.Options := [ofAllowMultiSelect,ofFileMustExist]

825

else

826

OpenHdrDlg.Options := [ofFileMustExist];

827

result := OpenHdrDlg.Execute;

828

if not result then exit;

829

if lForceMultiSelect then begin

830

lNumberofFiles:= OpenHdrDlg.Files.Count;

831

if lNumberofFiles < 2 then begin

832

Showmessage('Error: This function is designed to overlay MULTIPLE images. You selected less than two images.');

833

result := false;

834

end;

835

end;

836

end;

837

838

839

840

procedure TMainForm.NPMclick(Sender: TObject);

841

label

842

666;

843

var

844

lnGroup1,lMaskVoxels: integer;

845

lG: TStrings;

846

lMaskname: string;

847

lMaskHdr: TMRIcroHdr;

848

begin

849

if (not ttestmenu.checked) and (not BMmenu.checked) then begin

850

Showmessage('Error: you need to compute at least on test [options/test menu]');

851

exit;

852

end;

853

MsgClear;

854

Msg(GetKVers);

855

Msg('Threads: '+inttostr(gnCPUThreads));

856

if not OpenDialogExecute('Select brain mask ',false,false,kImgFilter) then begin

857

showmessage('NPM aborted: mask selection failed.');

858

exit;

859

end; //if not selected

860

lMaskname := OpenHdrDlg.Filename;

861

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

862

showmessage('Error reading mask.');

863

exit;

864

end;

865

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

866

if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file} then begin

867

showmessage('Mask file size too small.');

868

exit;

869

end;

870

Msg('Mask name = '+ lMaskname);

871

Msg('Total voxels = '+inttostr(lMaskVoxels));

872

//next, get 1st group

873

if not OpenDialogExecute('Select postive group (Z scores positive if this group is brighter)',true,true,kImgFilter) then begin

874

showmessage('NPM aborted: file selection failed.');

875

exit;

876

end; //if not selected

877

lG:= TStringList.Create; //not sure why TStrings.Create does not work???

878

lG.addstrings(OpenHdrDlg.Files);

879

lnGroup1 :=OpenHdrDlg.Files.Count;

880

Msg('Scans in Group 1 = '+inttostr(lnGroup1));

881

//next, get 2nd group

882

if not OpenDialogExecute('Select negative group (Z scores negative if this group is brighter)',true,true,kImgFilter) then begin

883

showmessage('NPM aborted: file selection failed.');

884

goto 666;

885

end; //if not selected

886

lG.addstrings(OpenHdrDlg.Files);

887

if not CheckVoxelsGroup(lG,lMaskVoxels) then begin

888

showmessage('File dimensions differ from mask.');

889

goto 666;

890

end;

891

Msg('Scans in Group 2 = '+inttostr(lG.count-lnGroup1));

892

NPMAnalyze(lG,lMaskHdr,lMaskVoxels,lnGroup1);

893

666:

894

lG.Free;

895

end;

896

897

function TMainForm.ThreshMap(lThresh: single; lVolVox: integer;lOutImg: singleP): integer;

898

var

899

lVox: integer;

900

begin

901

result := 0;

902

for lVox := 1 to lVolVox do

903

if lOutImg^[lVox] >= lThresh then

904

inc(result);

905

906

for lVox := 1 to lVolVox do

907

if lOutImg^[lVox] >= lThresh then

908

lOutImg^[lVox] := 1

909

else

910

lOutImg^[lVox] := 0;

911

end;

912

913

{x$DEFINE NOTmedianfx}

914

(*function TMainForm.LesionNPMAnalyze (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lnCrit,lRun: integer; var lSymptomRA: SingleP;var lFactname,lOutName: string): boolean;

915

label

916

123,667;

917

var

918

lOutNameMod: string;

919

lPlankImg: byteP;

920

lOutImgSum,lOutImgBM,lOutImgT,

921

lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP;

922

lPos,lPlank,lThread: integer;

923

lVolVox,lMinMask,lMaxMask,lTotalMemory,lnPlanks,lVoxPerPlank,

924

lThreadStart,lThreadEnd,lThreadInc,lnLesion,lnPermute,

925

lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lPosPct: int64;

926

lT,lBMz, lSum,lThresh :double;

927

lObsp: pointer;

928

lObs: Doublep0;

929

lStatHdr: TNIfTIhdr;

930

lFdata: file;

931

lRanOrderp: pointer;

932

lRanOrder: Doublep0;

933

lttest,lBM: boolean;

934

{$IFDEF medianfx}

935

lmedianFX,lmeanFX,lsummean,lsummedian: double;

936

lmediancount: integer;

937

{$ENDIF}

938

begin

939

lttest:= ttestmenu.checked;

940

lBM := BMmenu.checked;

941

lnPermute := ReadPermute;

942

Msg('Permutations = ' +IntToStr(lnPermute));

943

Msg('Analysis began = ' +TimeToStr(Now));

944

lTotalMemory := 0;

945

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

946

if (lVolVox < 1) then goto 667;

947

lMinMask := 1;

948

lMaxMask := lVolVox;

949

lVoxPerPlank := kPlankSz div lImages.Count div sizeof(byte) ;

950

if (lVoxPerPlank = 0) then goto 667; //no data

951

lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;

952

if (lTotalMemory = 0) then goto 667; //no data

953

lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;

954

Msg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));

955

Msg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));

956

if (lnPlanks = 1) then

957

getmem(lPlankImg,lTotalMemory) //assumes 1bpp

958

else

959

getmem(lPlankImg,kPlankSz);

960

lStartVox := lMinMask;

961

lEndVox := lMinMask-1;

962

{$IFDEF medianfx}

963

lsummean := 0;

964

lsummedian:= 0;

965

lmediancount := 0;

966

{$ENDIF}

967

for lPos := 1 to lImages.Count do

968

if gScaleRA[lPos] = 0 then

969

gScaleRA[lPos] := 1;

970

createArray64(lObsp,lObs,lImages.Count);

971

getmem(lOutImgSum,lVolVox* sizeof(single));

972

getmem(lOutImgBM,lVolVox* sizeof(single));

973

getmem(lOutImgT,lVolVox* sizeof(single));

974

InitPermute (lImages.Count, lnPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);

975

for lPos := 1 to lVolVox do begin

976

lOutImgSum^[lPos] := 0;

977

lOutImgBM^[lPos] := 0;

978

lOutImgT^[lPos] := 0;

979

end;

980

//next create permuted BM bounds

981

if lBM then begin

982

Msg('Generating BM permutation thresholds');

983

Refresh;

984

for lPos := 1 to lImages.Count do

985

lObs^[lPos-1] := lSymptomRA^[lPos];

986

genBMsim (lImages.Count, lObs);

987

end;

988

ClearThreadData(gnCPUThreads,lnPermute) ;

989

for lPlank := 1 to lnPlanks do begin

990

Msg('Computing plank = ' +Inttostr(lPlank));

991

Refresh;

992

Application.processmessages;

993

lEndVox := lEndVox + lVoxPerPlank;

994

if lEndVox > lMaxMask then begin

995

lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);

996

lEndVox := lMaxMask;

997

end;

998

lPlankImgPos := 1;

999

for lPos := 1 to lImages.Count do begin

1000

if not LoadImg8(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then

1001

goto 667;

1002

lPlankImgPos := lPlankImgPos + lVoxPerPlank;

1003

end;//for each image

1004

//threading start

1005

lThreadStart := 1;

1006

lThreadInc := lVoxPerPlank div gnCPUThreads;

1007

lThreadEnd := lThreadInc;

1008

Application.processmessages;

1009

for lThread := 1 to gnCPUThreads do begin

1010

if lThread = gnCPUThreads then

1011

lThreadEnd := lVoxPerPlank; //avoid integer rounding error

1012

with TLesionContinuous.Create (ProgressBar1,lttest,lBM,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,0,lPlankImg,lOutImgSum,lOutImgBM,lOutImgT,nil,lSymptomRA) do

1013

{$IFDEF FPC} OnTerminate := @ThreadDone; {$ELSE}OnTerminate := ThreadDone;{$ENDIF}

1014

inc(gThreadsRunning);

1015

lThreadStart := lThreadEnd + 1;

1016

lThreadEnd :=lThreadEnd + lThreadInc;

1017

end; //for each thread

1018

repeat

1019

Application.processmessages;

1020

until gThreadsRunning = 0;

1021

Application.processmessages;

1022

//threading end

1023

lStartVox := lEndVox + 1;

1024

end;

1025

lnVoxTested := SumThreadData(gnCPUThreads,lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);

1026

//next report findings

1027

if lnVoxTested < 1 then begin

1028

Msg('**Error: no voxels tested: no regions lesioned in at least '+inttostr(lnCrit)+' patients**');

1029

goto 123;

1030

end;

1031

1032

Msg('Voxels tested = ' +Inttostr(lnVoxTested));

1033

{$IFDEF medianfx}

1034

Msg('Average MEAN effect size = ' +realtostr((lsummean/lmediancount),3));

1035

Msg('Average MEDIAN effect size = ' +realtostr((lsummedian/lmediancount),3));

1036

{$ENDIF}

1037

Msg('Only tested voxels with more than '+inttostr(lnCrit)+' lesions');

1038

//Next: save results from permutation thresholding....

1039

reportBonferroni('Std',lnVoxTested);

1040

//next: save data

1041

MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);

1042

//save sum map

1043

lOutNameMod := ChangeFilePostfixExt(lOutName,'Sum'+lFactName,'.hdr');

1044

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);

1045

//create new header - subsequent images will use Z-scores

1046

MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );

1047

if Sum2PowerCont(lOutImgSum,lVolVox,lImages.Count) then begin

1048

lOutNameMod := ChangeFilePostfixExt(lOutName,'Power'+lFactName,'.hdr');

1049

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);

1050

end;

1051

1052

//MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );

1053

if lttest then begin //save Ttest

1054

//next: convert t-scores to z scores

1055

for lPos := 1 to lVolVox do

1056

lOutImgT^[lPos] := TtoZ (lOutImgT^[lPos],lImages.Count-2);

1057

for lPos := 1 to lnPermute do begin

1058

lPermuteMaxT^[lPos] := TtoZ (lPermuteMaxT^[lPos],lImages.Count-2);

1059

lPermuteMinT^[lPos] := TtoZ (lPermuteMinT^[lPos],lImages.Count-2);

1060

end;

1061

lThresh := reportFDR ('ttest', lVolVox, lnVoxTested, lOutImgT);

1062

reportPermute('ttest',lnPermute,lPermuteMaxT, lPermuteMinT);

1063

lOutNameMod := ChangeFilePostfixExt(lOutName,'ttest'+lFactName,'.hdr');

1064

if lRun > 0 then

1065

Msg('threshtt,'+inttostr(lRun)+','+inttostr(ThreshMap(lThresh,lVolVox,lOutImgT))+','+realtostr(lThresh,3));

1066

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgT,1);

1067

1068

end;

1069

if lBM then begin //save Mann Whitney

1070

lThresh := reportFDR ('BM', lVolVox, lnVoxTested, lOutImgBM);

1071

reportPermute('BM',lnPermute,lPermuteMaxBM, lPermuteMinBM);

1072

lOutNameMod := ChangeFilePostfixExt(lOutName,'BM'+lFactName,'.hdr');

1073

if lRun > 0 then

1074

Msg('threshbm,'+inttostr(lRun)+','+inttostr(ThreshMap(lThresh,lVolVox,lOutImgBM))+','+realtostr(lThresh,3));

1075

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgBM,1);

1076

1077

end;

1078

//next: free dynamic memory

1079

123:

1080

FreePermute (lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp);

1081

freemem(lOutImgT);

1082

freemem(lOutImgBM);

1083

freemem(lOutImgSum);

1084

freemem(lObsp);

1085

freemem(lPlankImg);

1086

Msg('Analysis finished = ' +TimeToStr(Now));

1087

lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes'+lFactName,'.txt');

1088

MsgSave(lOutNameMod);

1089

ProgressBar1.Position := 0;

1090

exit;

1091

667: //you only get here if you aborted ... free memory and report error

1092

if lTotalMemory > 1 then freemem(lPlankImg);

1093

Msg('Unable to complete analysis.');

1094

ProgressBar1.Position := 0;

1095

end; //LesionNPMAnalyze *)

1096

1097

1098

(*function TMainForm.LesionNPMAnalyzeBinomial (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lnCrit: integer; var lSymptomRA: SingleP; var lFactname,lOutName: string): boolean;

1099

label

1100

123,667;

1101

var

1102

lVal: single;

1103

lOutNameMod: string;

1104

lPlankImg: byteP;

1105

lOutImgSum,lOutImgL,lDummyImg,

1106

lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP;

1107

lPos,lPlank,lThread,lnDeficit: integer;

1108

lTotalMemory,lVolVox,lMinMask,lMaxMask,lnPlanks,lVoxPerPlank,

1109

lThreadStart,lThreadInc,lThreadEnd, lnLesion,lnPermute,

1110

lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lPosPct: int64;

1111

lT, lSum: double;

1112

lObsp: pointer;

1113

lObs: Doublep0;

1114

lStatHdr: TNIfTIhdr;

1115

lFdata: file;

1116

lRanOrderp: pointer;

1117

lRanOrder: Doublep0;

1118

begin

1119

lnPermute := ReadPermute;

1120

Msg('Permutations = ' +IntToStr(lnPermute));

1121

//lOutName := lMaskHdr.ImgFileName;

1122

//if not SaveHdrName ('Statistical Map', lOutName) then exit;

1123

Msg('Analysis began = ' +TimeToStr(Now));

1124

lTotalMemory := 0;

1125

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

1126

if (lVolVox < 1) then goto 667;

1127

lMinMask := 1;

1128

lMaxMask := lVolVox;

1129

lVoxPerPlank := kPlankSz div lImages.Count div sizeof(byte) ;

1130

if (lVoxPerPlank = 0) then goto 667; //no data

1131

lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;

1132

if (lTotalMemory = 0) then goto 667; //no data

1133

lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;

1134

Msg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));

1135

Msg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));

1136

if (lnPlanks = 1) then

1137

getmem(lPlankImg,lTotalMemory) //assumes 1bp

1138

else

1139

getmem(lPlankImg,kPlankSz);

1140

lStartVox := lMinMask;

1141

lEndVox := lMinMask-1;

1142

for lPos := 1 to lImages.Count do

1143

if gScaleRA[lPos] = 0 then

1144

gScaleRA[lPos] := 1;

1145

createArray64(lObsp,lObs,lImages.Count);

1146

getmem(lOutImgSum,lVolVox* sizeof(single));

1147

getmem(lOutImgL,lVolVox* sizeof(single));

1148

InitPermute (lImages.Count, lnPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);

1149

for lPos := 1 to lVolVox do begin

1150

lOutImgSum^[lPos] := 0;

1151

lOutImgL^[lPos] := 0;

1152

end;

1153

ClearThreadDataPvals(gnCPUThreads,lnPermute) ;

1154

for lPlank := 1 to lnPlanks do begin

1155

ProgressBar1.Position := 1;

1156

Msg('Computing plank = ' +Inttostr(lPlank));

1157

Refresh;

1158

Application.processmessages;

1159

lEndVox := lEndVox + lVoxPerPlank;

1160

if lEndVox > lMaxMask then begin

1161

lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);

1162

lEndVox := lMaxMask;

1163

end;

1164

lPlankImgPos := 1;

1165

for lPos := 1 to lImages.Count do begin

1166

if not LoadImg8(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then

1167

goto 667;

1168

lPlankImgPos := lPlankImgPos + lVoxPerPlank;

1169

end;//for each image

1170

//threading start

1171

lThreadStart := 1;

1172

lThreadInc := lVoxPerPlank div gnCPUThreads;

1173

lThreadEnd := lThreadInc;

1174

Application.processmessages;

1175

for lThread := 1 to gnCPUThreads do begin

1176

if lThread = gnCPUThreads then

1177

lThreadEnd := lVoxPerPlank; //avoid integer rounding error

1178

//with TLesionBinomial.Create (ProgressBar1,false,true,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,666, lDummyImg,lPlankImg,lOutImgSum,lOutImgL,lDummyImg,lSymptomRA) do

1179

with TLesionBinom.Create (ProgressBar1,false,true,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,0,lPlankImg,lOutImgSum,lOutImgL,lDummyImg,nil,lSymptomRA) do

1180

{$IFDEF FPC} OnTerminate := @ThreadDone; {$ELSE}OnTerminate := ThreadDone;{$ENDIF}

1181

inc(gThreadsRunning);

1182

Msg('Thread ' +Inttostr(gThreadsRunning)+' = '+inttostr(lThreadStart)+'..'+inttostr(lThreadEnd));

1183

lThreadStart := lThreadEnd + 1;

1184

lThreadEnd :=lThreadEnd + lThreadInc;

1185

end; //for each thread

1186

repeat

1187

Application.processmessages;

1188

until gThreadsRunning = 0;

1189

Application.processmessages;

1190

//threading end

1191

lStartVox := lEndVox + 1;

1192

end;

1193

lnVoxTested := SumThreadData(gnCPUThreads,lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);

1194

for lPos := 1 to lnPermute do begin

1195

if (lPermuteMinT^[lPos] > 1.1) or (lPermuteMinT^[lPos] < -1.1) then

1196

lPermuteMinT^[lPos] := 0.5;

1197

if (lPermuteMaxT^[lPos] > 1.1) or (lPermuteMaxT^[lPos] < -1.1) then

1198

lPermuteMaxT^[lPos] := 0.5;

1199

lVal := lPermuteMaxT^[lPos];

1200

lPermuteMaxT^[lPos] := lPermuteMinT^[lPos];

1201

lPermuteMinT^[lPos] := lVal;

1202

if lPermuteMaxT^[lPos] < 0 then

1203

lPermuteMaxT^[lPos] := -pNormalInv(abs(lPermuteMaxT^[lPos]))

1204

else

1205

lPermuteMaxT^[lPos] := pNormalInv(lPermuteMaxT^[lPos]);

1206

if lPermuteMinT^[lPos] < 0 then

1207

lPermuteMinT^[lPos] := -pNormalInv(abs(lPermuteMinT^[lPos]))

1208

else

1209

lPermuteMinT^[lPos] := pNormalInv(lPermuteMinT^[lPos]);

1210

end;

1211

1212

1213

1214

if lnVoxTested < 1 then begin

1215

Msg('**Error: no voxels tested: no regions lesioned in at least '+inttostr(lnCrit)+' patients**');

1216

goto 123;

1217

end;

1218

//next report findings

1219

Msg('Voxels tested = ' +Inttostr(lnVoxTested));

1220

Msg('Only tested voxels with more than '+inttostr(lnCrit)+' lesions');

1221

//Next: save results from permutation thresholding....

1222

reportBonferroni('Std',lnVoxTested);

1223

//next: save data

1224

//savedata

1225

MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);

1226

//save sum map

1227

lOutNameMod := ChangeFilePostfixExt(lOutName,'Sum'+lFactName,'.hdr');

1228

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);

1229

//future images will store Z-scores...

1230

MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );

1231

//save power map

1232

lnDeficit := 0;

1233

for lPos := 1 to lImages.Count do

1234

if lSymptomRA^[lPos] = 0 then

1235

inc(lnDeficit);

1236

if Sum2Power(lOutImgSum,lVolVox,lImages.Count,lnDeficit) then begin

1237

lOutNameMod := ChangeFilePostfixExt(lOutName,'Power'+lFactName,'.hdr');

1238

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);

1239

end;

1240

// MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );

1241

//save Liebermeister

1242

1243

lOutNameMod := ChangeFilePostfixExt(lOutName,'L'+lFactName,'.hdr');

1244

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgL,1);

1245

//save end

1246

reportFDR ('L', lVolVox, lnVoxTested, lOutImgL);

1247

reportPermute('L',lnPermute,lPermuteMaxT, lPermuteMinT);

1248

1249

123:

1250

//next: free dynamic memory

1251

FreePermute (lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp);

1252

freemem(lOutImgL);

1253

freemem(lOutImgSum);

1254

freemem(lObsp);

1255

freemem(lPlankImg);

1256

Msg('Analysis finished = ' +TimeToStr(Now));

1257

lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes'+lFactName,'.txt');

1258

MsgSave(lOutNameMod);

1259

1260

ProgressBar1.Position := 0;

1261

exit;

1262

667: //you only get here if you aborted ... free memory and report error

1263

if lTotalMemory > 1 then freemem(lPlankImg);

1264

Msg('Unable to complete analysis.');

1265

ProgressBar1.Position := 0;

1266

end; *)

1267

1268

1269

function TMainForm.GetValX (var lnSubj, lnFactors: integer; var lSymptomRA: singleP; var lImageNames: TStrings; var lCrit: integer; {lBinomial : boolean;} var lPredictorList: TStringList):boolean;

1270

//warning: you MUST free lPredictorList

1271

var

1272

lVALFilename {,lTemplateName}: string;

1273

lCritPct: integer;

1274

begin

1275

lPredictorList := TStringList.Create;

1276

result := false;

1277

lnSubj := 0;

1278

if not MainForm.OpenDialogExecute('Select MRIcron VAL file',false,false,'MRIcron VAL (*.val)|*.val') then begin

1279

showmessage('NPM aborted: VAL file selection failed.');

1280

exit;

1281

end; //if not selected

1282

lVALFilename := MainForm.OpenHdrDlg.Filename;

1283

result := GetValCore ( lVALFilename, lnSubj, lnFactors, lSymptomRA, lImageNames, lCrit,lCritPct{,binom},lPredictorList);

1284

end;

1285

1286

1287

function TMainForm.ReportDescriptives (var RA: SingleP; n: integer): boolean;

1288

var lMn,lSD,lSE,lSkew,lZSkew: double;

1289

begin

1290

SuperDescriptive (RA, n, lMn,lSD,lSE,lSkew,lZSkew);

1291

Msg('mean='+floattostr(lMn)+',StDev='+floattostr(lSD)+',StEr='+floattostr(lSE)+',Skew='+floattostr(lSkew)+',ZSkew='+floattostr(lZSkew));

1292

end;

1293

1294

(*function noVariance (lRA: singlep; lnSubj: integer): boolean;

1295

var

1296

lI : integer;

1297

begin

1298

result := false;

1299

if lnSubj < 2 then exit;

1300

for lI := 2 to lnSubj do

1301

if lRA^[1] <> lRA^[lI] then

1302

exit;

1303

result := true;

1304

end; *)

1305

1306

(*procedure TMainForm.LesionBtnClick(Sender: TObject);

1307

label

1308

666;

1309

var

1310

lBinomial: boolean;

1311

lFact,lnFactors,lSubj,lnSubj,lnSubjAll,lMaskVoxels,lnCrit: integer;

1312

lImageNames,lImageNamesAll: TStrings;

1313

lPredictorList: TStringList;

1314

lTemp4D,lMaskname,lOutName,lFactname: string;

1315

lMaskHdr: TMRIcroHdr;

1316

lMultiSymptomRA,lSymptomRA: singleP;

1317

begin

1318

lBinomial := not odd( (Sender as tMenuItem).tag);

1319

if (not lBinomial) and (not ttestmenu.checked) and (not BMmenu.checked) then begin

1320

Showmessage('Error: you need to compute at least on test [options/test menu]');

1321

exit;

1322

end;

1323

lImageNamesAll:= TStringList.Create; //not sure why TStrings.Create does not work???

1324

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

1325

//next, get 1st group

1326

if not GetVal(lnSubjAll,lnFactors,lMultiSymptomRA,lImageNamesAll,lnCrit{,binom},lPredictorList) then begin

1327

showmessage('Error with VAL file');

1328

goto 666;

1329

end;

1330

lTemp4D := CreateDecompressed4D(lImageNamesAll);

1331

if (lnSubjAll < 1) or (lnFactors < 1) then begin

1332

Showmessage('Not enough subjects ('+inttostr(lnSubjAll)+') or factors ('+inttostr(lnFactors)+').');

1333

goto 666;

1334

end;

1335

lMaskname := lImageNamesAll[0];

1336

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

1337

showmessage('Error reading 1st mask.');

1338

goto 666;

1339

end;

1340

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

1341

if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file} then begin

1342

showmessage('Mask file size too small.');

1343

goto 666;

1344

end;

1345

if not CheckVoxelsGroup(lImageNamesAll,lMaskVoxels) then begin

1346

showmessage('File dimensions differ from mask.');

1347

goto 666;

1348

end;

1349

lOutName := ExtractFileDirWithPathDelim(lMaskName)+'results';

1350

SaveHdrDlg.Filename := loutname;

1351

lOutName := lOutName+'.nii.gz';

1352

if not SaveHdrName ('Base Statistical Map', lOutName) then goto 666;

1353

for lFact := 1 to lnFactors do begin

1354

MsgClear;

1355

Msg(GetKVers);

1356

lImageNames.clear;

1357

for lSubj := 1 to lnSubjAll do

1358

if (not lBinomial) or (lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)] = 0) OR (lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)] = 1) THEN begin

1359

{$IFNDEF FPC}if lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)] <> NaN then {$ENDIF}

1360

lImageNames.Add(lImageNamesAll[lSubj-1]);

1361

end else begin

1362

Msg('Data rejected: behavior must be zero or one for binomial test '+lImageNamesAll.Strings[lSubj-1]);

1363

end;

1364

lnSubj := lImageNames.Count;

1365

if lnSubj > 1 then begin

1366

getmem(lSymptomRA,lnSubj * sizeof(single));

1367

lnSubj := 0;

1368

for lSubj := 1 to lnSubjAll do

1369

if (not lBinomial) or (lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)] = 0) OR (lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)] = 1) THEN

1370

{$IFNDEF FPC}if lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)] <> NaN then begin

1371

{$ELSE} begin{$ENDIF}

1372

inc(lnSubj);

1373

lSymptomRA^[lnSubj] := lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)];

1374

end;

1375

Msg('Threads: '+inttostr(gnCPUThreads));

1376

lFactName := lPredictorList.Strings[lFact-1];

1377

lFactName := LegitFilename(lFactName,lFact);

1378

Msg('Factor = '+lFactname);

1379

For lSubj := 1 to lnSubj do

1380

Msg (lImageNames.Strings[lSubj-1] + ' = '+realtostr(lSymptomRA^[lSubj],2) );

1381

Msg('Total voxels = '+inttostr(lMaskVoxels));

1382

Msg('Only testing voxels damaged in at least '+inttostr(lnCrit)+' individual[s]');

1383

Msg('Number of Lesion maps = '+inttostr(lnSubj));

1384

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

1385

showmessage('File dimensions differ from mask.');

1386

goto 666;

1387

end;

1388

if noVariance (lSymptomRA,lnSubj) then

1389

Msg('Error no variability in behavioral data ')

1390

else if lBinomial then

1391

LesionNPMAnalyzeBinomial2(lImageNames,lMaskHdr,lnCrit,MainForm.ReadPermute,lSymptomRA,lFactname,lOutName)

1392

else begin

1393

ReportDescriptives(lSymptomRA,lnSubj);

1394

LesionNPMAnalyze2(lImageNames,lMaskHdr,lnCrit,-1,MainForm.ReadPermute,lSymptomRA,lFactName,lOutname,ttestmenu.checked,BMmenu.checked);

1395

end;

1396

Freemem(lSymptomRA);

1397

end; //lnsubj > 1

1398

end; //for each factor

1399

if lnSubjAll > 0 then begin

1400

Freemem(lMultiSymptomRA);

1401

end;

1402

666:

1403

lImageNames.Free;

1404

lImageNamesAll.Free;

1405

lPredictorList.Free;

1406

DeleteDecompressed4D(lTemp4D);

1407

end; *)

1408

1409

procedure TMainForm.Copy1Click(Sender: TObject);

1410

begin

1411

Memo1.SelectAll;

1412

Memo1.CopyToClipBoard;

1413

end;

1414

1415

procedure TMainForm.testmenuclick(Sender: TObject);

1416

begin

1417

(sender as TMenuItem).checked := not (sender as TMenuItem).checked;

1418

end;

1419

1420

procedure TMainForm.radiomenuclick(Sender: TObject);

1421

begin

1422

(sender as tmenuitem).checked := true;

1423

end;

1424

1425

procedure ComputePlankSize;

1426

begin

1427

if kPlankMB < 128 then

1428

kPlankMB := 128;

1429

if kPlankMB > 2000 then

1430

kPlankMB := 2000; //we use signed 32-bit pointers, so we can not exceed 2Gb

1431

kPlankSz :=1024 {bytes/kb} * 1024 {bytes/mb} * kPlankMB;

1432

kVers := kVers + ' CacheMB = '+inttostr(kPlankMB);

1433

end;

1434

1435

procedure TMainForm.ReadIniFile;

1436

var

1437

lFilename: string;

1438

lThreads: integer;

1439

lIniFile: TIniFile;

1440

begin

1441

lFilename := IniName;

1442

if not FileexistsEx(lFilename) then

1443

exit;

1444

lIniFile := TIniFile.Create(lFilename);

1445

ttestmenu.checked := IniBool(lIniFile,'computettest',true);

1446

//welchmenu.checked := IniBool(lIniFile,'computewelch',true);

1447

BMmenu.checked := IniBool(lIniFile,'computebm',false);

1448

gNULP := IniBool(lIniFile,'countlesionpatterns',false);

1449

gROI := IniBool(lIniFile,'ROI',false);

1450

kPlankMB := IniInt(lIniFile,'CacheMB',512);

1451

1452

WritePermute(IniInt(lIniFile,'nPermute',0));

1453

lThreads := IniInt(lIniFile,'nThread', gnCPUThreads );

1454

if lThreads > gnCPUThreads then

1455

lThreads := gnCPUThreads;

1456

gnCPUThreads := lThreads;

1457

lIniFile.Free;

1458

end; //ReadIniFile

1459

1460

procedure TMainForm.WriteIniFile;

1461

var

1462

lIniName: string;

1463

lIniFile: TIniFile;

1464

begin

1465

//showmessage('aaa');

1466

lIniName := IniName;

1467

if (DiskFreeEx(lIniName) < 1) then

1468

exit;

1469

lIniFile := TIniFile.Create(lIniName);

1470

lIniFile.WriteString('BOOL', 'computettest',Bool2Char(ttestmenu.checked));

1471

lIniFile.WriteString('BOOL', 'countlesionpatterns',Bool2Char(gNULP));

1472

lIniFile.WriteString('BOOL', 'ROI',Bool2Char(gROI));

1473

1474

//lIniFile.WriteString('BOOL', 'computewelch',Bool2Char(welchmenu.checked));

1475

lIniFile.WriteString('BOOL', 'computebm',Bool2Char(BMmenu.checked));

1476

lIniFile.WriteString('INT', 'CacheMB',inttostr(kPlankMB));

1477

lIniFile.WriteString('INT', 'nPermute',inttostr(ReadPermute));

1478

lIniFile.WriteString('INT', 'nThread',inttostr(ReadThread));

1479

lIniFile.Free;

1480

end;

1481

1482

1483

1484

procedure TMainForm.FormCreate(Sender: TObject);

1485

begin

1486

{$IFDEF Darwin}

1487

File1.visible := false;//for OSX, exit is in the application's menu

1488

Edit1.visible := false;//clipboard note yet working for OSX

1489

{$ENDIF}

1490

gnCPUThreads := GetLogicalCpuCount;

1491

(*if (ssShift in KeyDataToShiftState(vk_Shift)) then begin

1492

case MessageDlg('Shift key down during launch: do you want to reset the default preferences?', mtConfirmation,

1493

[mbYes, mbNo], 0) of { produce the message dialog box }

1494

mrNo: ReadIniFile;

1495

end; //case

1496

end else *)

1497

if not ResetDefaults then

1498

ReadIniFile;

1499

WriteThread(gnCPUThreads);

1500

ComputePlankSize;

1501

ROIanalysis1.visible := gROI;

1502

{$IFDEF compileANACOM}

1503

AnaCOMmenu.visible := gROI;

1504

{$ENDIF}

1505

end;

1506

1507

1508

1509

1510

function TMainForm.MakeMean (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lBinarize,lVariance : boolean): boolean;

1511

label

1512

667;

1513

var

1514

lOutName,lOutNameMod: string;

1515

lCountRA,lOutImgMn,lOutStDev,lPlankImg: SingleP;

1516

lTotalMemory: double;

1517

lPlank,lVolVox,lPos,lMinMask,lMaxMask,lnPlanks,lVoxPerPlank,

1518

lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lPosPct: integer;

1519

lStDev: boolean;

1520

lT, lSum,lSumSqr,lSD, lMn,lTotalSum,lTotalN: double;

1521

lStatHdr: TNIfTIhdr;

1522

lFdata: file;

1523

begin

1524

result := false;

1525

if not SaveHdrName ('Output image', lOutName) then exit;

1526

if (not lVariance) and (not lBinarize) then

1527

lStDev := true

1528

else

1529

lStDev := false;

1530

if lStDev then

1531

lStDev := OKMsg('Create a standard deviation image as well as a mean image?');

1532

Msg('Analysis began = ' +TimeToStr(Now));

1533

lTotalMemory := 0;

1534

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

1535

if (lVolVox < 1) then goto 667;

1536

lMinMask := 1;

1537

lMaxMask := lVolVox;

1538

lVoxPerPlank := kPlankSz div lImages.Count div sizeof(single) ;

1539

if (lVoxPerPlank = 0) then goto 667; //no data

1540

lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;

1541

if (lTotalMemory = 0) then goto 667; //no data

1542

lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;

1543

Msg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));

1544

Msg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));

1545

// fx(kPlankSz,8888);

1546

getmem(lPlankImg,kPlankSz);

1547

lStartVox := lMinMask;

1548

lEndVox := lMinMask-1;

1549

Msg('Number of scans = '+inttostr(lImages.count));

1550

Msg(' Index,Filename,Intercept,Slope');

1551

if lBinarize then begin

1552

getmem(lCountRA,lImages.Count*sizeof(single));

1553

for lPos := 1 to lImages.Count do begin

1554

gInterceptRA[lPos] := 0;

1555

gScaleRA[lPos] := 1;

1556

lCountRA^[lPos] := 0;

1557

end;

1558

end else begin

1559

for lPos := 1 to lImages.Count do begin

1560

Msg(' '+inttostr(lPos)+','+lImages[lPos-1]+','+realtostr(gInterceptRA[lPos],4)+','+realtostr(gScaleRA[lPos],4));

1561

if gScaleRA[lPos] = 0 then

1562

gScaleRA[lPos] := 1;

1563

end;

1564

end;

1565

1566

lTotalSum := 0;

1567

lTotalN := 0;

1568

//createArray64(lObsp,lObs,lImages.Count);

1569

getmem(lOutImgMn,lVolVox* sizeof(single));

1570

if lStDev then

1571

getmem(lOutStDev,lVolVox* sizeof(single));

1572

for lPlank := 1 to lnPlanks do begin

1573

Msg('Computing plank = ' +Inttostr(lPlank));

1574

Refresh;

1575

lEndVox := lEndVox + lVoxPerPlank;

1576

if lEndVox > lMaxMask then begin

1577

lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);

1578

lEndVox := lMaxMask;

1579

end;

1580

lPlankImgPos := 1;

1581

for lPos := 1 to lImages.Count do begin

1582

if not LoadImg(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then

1583

goto 667;

1584

lPlankImgPos := lPlankImgPos + lVoxPerPlank;

1585

end;//for each image

1586

lPosPct := lVoxPerPlank div 100;

1587

for lPos2 := 1 to lVoxPerPlank do begin

1588

if (lPos2 mod lPosPct) = 0 then begin

1589

ProgressBar1.Position := round((lPos2/lVoxPerPlank)*100);

1590

Application.Processmessages;

1591

end;

1592

lPos2Offset := lPos2+lStartVox-1;

1593

lSum := 0;

1594

if lVariance then begin

1595

lSum := sqr(lPlankImg^[lPos2]-lPlankImg^[lVoxPerPlank+lPos2]);//actually variance...

1596

//% signal

1597

//if lPlankImg[lVoxPerPlank+lPos2] <> 0 then

1598

1599

// lSum := lPlankImg[lPos2]/lPlankImg[lVoxPerPlank+lPos2]

1600

//else

1601

// lSum := 0;//pct signal...

1602

//end % signal

1603

lOutImgMn^[lPos2Offset] := lSum;

1604

lTotalSum := lTotalSum + lOutImgMn^[lPos2Offset];

1605

lTotalN := lTotalN + 1;

1606

end else begin //not variance

1607

1608

if lBinarize then begin

1609

for lPos := 1 to lImages.Count do

1610

if lPlankImg^[((lPos-1)* lVoxPerPlank)+lPos2] <> 0 then begin

1611

lSum := lSum+1;

1612

lCountRA^[lPos] := lCountRA^[lPos] + 1;

1613

end;

1614

end else

1615

for lPos := 1 to lImages.Count do

1616

lSum := lSum +(gScaleRA[lPos]*lPlankImg^[((lPos-1)* lVoxPerPlank)+lPos2])+gInterceptRA[lPos];

1617

// fx(lPos, gScaleRA[lPos],gInterceptRA[lPos]);

1618

lOutImgMn^[lPos2Offset] := lSum/lImages.Count;

1619

if lStDev then begin

1620

//lSum := 0;

1621

//for lPos := 1 to lImages.Count do

1622

// lSum := lSum + (gScaleRA[lPos]*lPlankImg^[((lPos-1)* lVoxPerPlank)+lPos2])+gInterceptRA[lPos];

1623

lSumSqr := 0;

1624

for lPos := 1 to lImages.Count do

1625

lSumSqr := lSumSqr + Sqr((gScaleRA[lPos]*lPlankImg^[((lPos-1)* lVoxPerPlank)+lPos2])+gInterceptRA[lPos]);

1626

lSD := (lSumSqr - ((Sqr(lSum))/lImages.Count));

1627

if (lSD > 0) then

1628

lSD := Sqrt ( lSD/(lImages.Count-1))

1629

else begin

1630

lSD := 0;

1631

(*if l1stError then begin

1632

for lPos := 1 to lImages.Count do

1633

Msg(floattostr( (gScaleRA[lPos]*lPlankImg^[((lPos-1)* lVoxPerPlank)+lPos2])+gInterceptRA[lPos]));

1634

msg('---');

1635

msg(floattostr(lSum));

1636

msg(floattostr(lSumSqr));

1637

1638

l1stError := false;

1639

end;*)

1640

end;

1641

lOutStDev^[lPos2Offset] := lSD;

1642

end;

1643

end; //not variance

1644

if lSum > 0 then begin

1645

lTotalSum := lTotalSum + lOutImgMn^[lPos2Offset];

1646

lTotalN := lTotalN + 1;

1647

end;

1648

1649

end;

1650

lStartVox := lEndVox + 1;

1651

end;

1652

if lBinarize then begin

1653

for lPos := 1 to lImages.Count do begin

1654

Msg(' '+inttostr(lPos)+','+lImages[lPos-1]+','+inttostr(round(lCountRA^[lPos])) );

1655

1656

lCountRA^[lPos] := 0;

1657

end;

1658

freemem(lCountRA);

1659

end; //if binar

1660

//next: save data

1661

MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);

1662

//save mean

1663

1664

1665

if lVariance then

1666

lOutNameMod := ChangeFilePostfixExt(lOutName,'var','.hdr')

1667

else

1668

lOutNameMod := ChangeFilePostfixExt(lOutName,'Mean','.hdr');

1669

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgMn,1);

1670

freemem(lOutImgMn);

1671

if lStDev then begin

1672

lOutNameMod := ChangeFilePostfixExt(lOutName,'StDev','.hdr');

1673

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutStDev,1);

1674

freemem(lOutStDev);

1675

end;

1676

1677

//freemem(lObsp);

1678

freemem(lPlankImg);

1679

Msg('Analysis finished = ' +TimeToStr(Now));

1680

lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes','.txt');

1681

MsgSave(lOutNameMod);

1682

if (lTotalN > 0) then

1683

Msg('num voxels >0 = ' +inttostr(round(lTotalN))+' mean value for voxels >0: '+floattostr(lTotalSum/lTotalN));

1684

1685

ProgressBar1.Position := 0;

1686

exit;

1687

667: //you only get here if you aborted ... free memory and report error

1688

if lTotalMemory > 1 then freemem(lPlankImg);

1689

Msg('Unable to complete analysis.');

1690

ProgressBar1.Position := 0;

1691

end;

1692

1693

procedure TMainForm.Makemeanimage1Click(Sender: TObject);

1694

label

1695

666;

1696

var

1697

lMaskVoxels: integer;

1698

lG: TStrings;

1699

lMaskname: string;

1700

lMaskHdr: TMRIcroHdr;

1701

begin

1702

MsgClear;

1703

Msg(GetKVers);

1704

if not OpenDialogExecute('Select images to average',true,true,kImgFilter) then begin

1705

showmessage('NPM aborted: file selection failed.');

1706

exit;

1707

end; //if not selected

1708

lG:= TStringList.Create;

1709

lG.addstrings(OpenHdrDlg.Files);

1710

lMaskname := lG[0];

1711

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

1712

showmessage('Error reading '+lMaskName);

1713

goto 666;

1714

end;

1715

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

1716

if not CheckVoxelsGroup(lG,lMaskVoxels) then begin

1717

showmessage('File dimensions differ from mask.');

1718

goto 666;

1719

end;

1720

Msg('Voxels = '+inttostr(lMaskVoxels));

1721

MakeMean(lG,lMaskHdr, odd((Sender as TMenuItem).tag),false);

1722

666:

1723

lG.Free;

1724

end;

1725

1726

procedure TMainForm.Exit1Click(Sender: TObject);

1727

begin

1728

Close;

1729

end;

1730

1731

function MinMax (var lImg: SingleP; var lVolVox: integer; var lMin, lMax: single): boolean;

1732

var

1733

lC: integer;

1734

begin

1735

result := false;

1736

if lVolVox < 1 then

1737

exit;

1738

lMax := lImg^[1];

1739

for lC := 1 to lVolVox do

1740

if lImg^[lC] > lMax then

1741

lMax := lImg^[lC];

1742

//lCx := lC;

1743

lMin := lImg^[1];

1744

for lC := 1 to lVolVox do

1745

if lImg^[lC] < lMin then

1746

lMin := lImg^[lC];

1747

result := true;

1748

end;

1749

1750

function DetectMode (var lImg: SingleP; var lVolVox: integer; var lMin, lMax, lModeLo,lModeHi: single; lInflection: boolean): boolean;

1751

const

1752

kHistoBins = 255;//numbers of bins for histogram/image balance

1753

var

1754

lSmooth,lPrevSmooth,lModeWid,lC,lMinPos,lMode,lModePos,lMaxModePos,lMode2NotInflection: integer;

1755

lMod,lRng: single;

1756

lHisto : array [0..kHistoBins] of longint;

1757

begin

1758

1759

result := false;

1760

if (lVolVox < 1) or (lMax < lMin) then

1761

exit;

1762

//zero array

1763

for lC := 1 to kHistoBins do

1764

lHisto[lC] := 0;

1765

//find scaling

1766

lRng := abs(lMax-lMin);

1767

if lRng > 0 then

1768

lMod := (kHistoBins)/lRng

1769

else

1770

lMod := 0;

1771

//fill histogram

1772

for lC := 1 to lVolVox do

1773

if lImg^[lC] <> 0 then

1774

inc(lHisto[round((lImg^[lC]-lMin)*lMod)]);

1775

1776

{for lC := 1 to lVolVox do

1777

inc(lHisto[round((lImg^[lC]-lMin)*lMod)]); }

1778

//smooth histogram

1779

lPrevSmooth := lHisto[1];

1780

for lC := 2 to (kHistoBins-1) do begin

1781

lSmooth := round( (lHisto[lC-1]+lHisto[lC]+lHisto[lC]+lHisto[lC+1])/4);

1782

lHisto[lC-1] := lPrevSmooth;

1783

lPrevSmooth := lSmooth;

1784

end;

1785

lHisto[kHistoBins-1] := lPrevSmooth;

1786

//find mode

1787

lMode := 0;

1788

lMinPos := 1;//indexed from zero

1789

//find highest peak

1790

for lC := lMinPos to kHistoBins do begin

1791

if lHisto[lC] > lMode then begin

1792

lModePos := lC;

1793

lMode := lHisto[lC];

1794

end;//if new mode

1795

end; //for each bin

1796

if lMode > 0 then

1797

lMaxModePos := lModePos

1798

else

1799

exit;

1800

//find 2nd highest peak

1801

//find 2nd highest peak

1802

lModeWid := 25;

1803

lModePos := lMinPos;

1804

lMode := lHisto[lMinPos];

1805

if (lMaxModePos - lModeWid) > lMinPos then begin

1806

for lC := lMinPos to (lMaxModePos - lModeWid) do begin

1807

if lHisto[lC] > lMode then begin

1808

lModePos := lC;

1809

lMode := lHisto[lC];

1810

end;//if new mode

1811

end; //for each bin

1812

end; //check below highest peak

1813

if (lMaxModePos + lModeWid) < kHistoBins then begin

1814

for lC := (lMaxModePos + lModeWid) to kHistoBins do begin

1815

if lHisto[lC] > lMode then begin

1816

lModePos := lC;

1817

lMode := lHisto[lC];

1818

end;//if new mode

1819

end; //for each bin

1820

end; //check above highest peak

1821

//fx(lModePos);

1822

//an alternative method to find mode is to look for inflection - less assumptions, more sensitive to noise

1823

if lInflection then begin

1824

lMode2NotInflection := lModePos;

1825

lModePos := lMinPos;

1826

1827

lMode := 0;

1828

lC := lMaxModePos;

1829

while ((lC-1) > lMinPos) and (lHisto[lC] > lHisto[lC-1]) do

1830

dec(lC); //find inflection

1831

while ((lC-1) > lMinPos) do begin

1832

dec(lC);

1833

if lHisto[lC] > lMode then begin

1834

lModePos := lC;

1835

lMode := lHisto[lC];

1836

end;//if new mode

1837

end; //look for mode

1838

1839

lC := lMaxModePos;

1840

while ((lC+1) <= kHistoBins) and (lHisto[lC] > lHisto[lC+1]) do

1841

inc(lC); //find inflection

1842

while ((lC+1) <= kHistoBins) do begin

1843

inc(lC);

1844

if lHisto[lC] > lMode then begin

1845

lModePos := lC;

1846

lMode := lHisto[lC];

1847

end;//if new mode

1848

end; //look for mode

1849

1850

if abs(lMode2NotInflection-lModePos) > 3 then

1851

Showmessage('Warning: inflection and windowed algorithms find different 2nd modes. Using inflection 2nd mode. inflection ='+inttostr(lModePos)+' windowed: '+inttostr(lMode2NotInflection));

1852

1853

end;

1854

//now, return scaled values...

1855

if lMod = 0 then exit;

1856

lModeLo := (lModePos/lMod)+lMin;

1857

lModeHi := (lMaxModePos/lMod)+lMin;

1858

if lModeLo > lModeHi then begin

1859

lMod := lModeLo;

1860

lModeLo := lModeHi;

1861

lModeHi := lMod;

1862

end;

1863

result := true;

1864

end;

1865

1866

1867

procedure CopyFileEXoverwrite (lInName,lOutName: string);

1868

var lFSize: Integer;

1869

lBuff: bytep0;

1870

lFData: file;

1871

begin

1872

lFSize := FSize(lInName);

1873

if (lFSize < 1) then exit;

1874

assignfile(lFdata,lInName);

1875

filemode := 0;

1876

reset(lFdata,lFSize{1});

1877

GetMem( lBuff, lFSize);

1878

BlockRead(lFdata, lBuff^, 1{lFSize});

1879

closefile(lFdata);

1880

assignfile(lFdata,lOutName);

1881

filemode := 2;

1882

Rewrite(lFdata,lFSize);

1883

BlockWrite(lFdata,lBuff^, 1 {, NumWritten});

1884

closefile(lFdata);

1885

freemem(lBuff);

1886

end;

1887

1888

procedure CopyFileEX (lInName,lOutName: string);

1889

var lFSize: Integer;

1890

begin

1891

lFSize := FSize(lInName);

1892

if (lFSize < 1) or (fileexistsEX(lOutName)) then exit;

1893

CopyFileEXoverwrite (lInName,lOutName);

1894

end;

1895

1896

1897

function DetectMeanStDev (var lImg: SingleP; var lVolVox: integer; var lMean,lStDev: double): boolean;

1898

var

1899

lIncVox,lVox: integer;

1900

lSum,lSumSqr,lSE: double;

1901

begin

1902

lMean := 0;

1903

lStDev := 0;

1904

result := false;

1905

if (lVolVox < 1) then

1906

exit;

1907

lSum := 0;

1908

lSumSqr := 0;

1909

lIncVox := 0; //voxels included - e.g. not masked

1910

for lVox := 1 to lVolVox do begin

1911

if lImg^[lVox] <> 0 then begin //not in mask

1912

inc(lIncVox);

1913

lSum := lSum + lImg^[lVox];

1914

lSumSqr := lSumSqr + sqr(lImg^[lVox]);

1915

end;

1916

end;

1917

if lincVox < 3 then

1918

exit;

1919

Descriptive (lincVox, lSumSqr, lSum,lMean,lStDev,lSE);

1920

result := true;

1921

end; //DetectMeanStDev

1922

//zero array

1923

1924

1925

1926

function TMainForm.Balance (var lImageName,lMaskName: String; {lInflection: boolean}lMethod: integer): boolean;

1927

//0 = masked peak

1928

//1 = inflection

1929

//2 = mean =1, stdev=1

1930

var

1931

lImg,lMaskImg: SingleP;

1932

lHdr,lMaskHdr: TMRIcroHdr;

1933

lVolVox,lVox,lMasked: integer;

1934

lMaskedInten,lMean,lStDev: double;

1935

lMin,lMax: single;

1936

lModeLo,lModeHi,lIntercept,lSlope: single;

1937

lOutNameMod: string;

1938

begin

1939

//lOutName := lMaskHdr.ImgFileName;

1940

result := false;

1941

//if not SaveHdrName ('Statistical Map', lOutName) then exit;

1942

if not NIFTIhdr_LoadHdr(lImageName,lHdr) then begin

1943

showmessage('Error reading '+lImageName);

1944

exit;

1945

end;

1946

lVolVox := lHdr.NIFTIhdr.dim[1]*lHdr.NIFTIhdr.dim[2]* lHdr.NIFTIhdr.dim[3];

1947

if (lVolVox < 1) then exit;

1948

getmem(lImg,lVolVox*sizeof(single));

1949

if not LoadImg(lHdr.ImgFileName, lImg, 1, lVolVox,round(lHdr.NIFTIhdr.vox_offset),1,lHdr.NIFTIhdr.datatype,lVolVox) then begin

1950

Msg('Unable to load ' +lHdr.ImgFileName);

1951

exit;

1952

end;

1953

if lMaskName <> '' then begin

1954

if not NIFTIhdr_LoadHdr(lMaskName,lMaskHdr) then begin

1955

showmessage('Error reading '+lMaskName);

1956

exit;

1957

end;

1958

if lVolVox <> (lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3]) then begin

1959

showmessage('Mask and header must have identical dimensions '+lMaskName+ ' ' + lImageName);

1960

exit;

1961

1962

end;

1963

getmem(lMaskImg,lVolVox*sizeof(single));

1964

if not LoadImg(lMaskHdr.ImgFileName, lMaskImg, 1, lVolVox,round(lMaskHdr.NIFTIhdr.vox_offset),1,lMaskHdr.NIFTIhdr.datatype,lVolVox) then begin

1965

Msg('Unable to load mask ' +lMaskHdr.ImgFileName);

1966

exit;

1967

end;

1968

lMasked := 0;

1969

lMaskedInten := 0;

1970

for lVox := 1 to lVolVox do

1971

if lMaskImg^[lVox] = 0 then begin

1972

lMaskedInten := lMaskedInten + lImg^[lVox];

1973

lImg^[lVox] := 0;

1974

inc(lMasked);

1975

end;

1976

if lMasked < 1 then

1977

Msg('Warning: no voxels masked with image '+lMaskName)

1978

else

1979

Msg('Mask='+ lMaskName+' Number of voxels masked= '+inttostr(lMasked)+' Mean unscaled intensity of masked voxels= '+floattostr(lMaskedInten/lMasked));

1980

freemem(lMaskImg);

1981

end;//mask

1982

1983

if not MinMax(lImg,lVolVox,lMin,lMax) then exit;

1984

Msg(lImageName+' -> '+lHdr.ImgFileName);

1985

Msg('min..max ' +floattostr(lMin)+'..'+floattostr(lMax));

1986

if (lMethod = 0) or (lMethod = 1) then begin

1987

if not DetectMode(lImg,lVolVox,lMin,lMax,lModeLo,lModeHi, odd(lMethod)) then exit;

1988

if odd(lMethod) then

1989

Msg('method for finding second mode: inflection')

1990

else

1991

Msg('method for finding second mode: masked peak');

1992

Msg('modes Lo Hi ' +floattostr(lModeLo)+'..'+floattostr(lModeHi));

1993

if lModeLo >= lModeHi then exit;

1994

lSlope := 1/abs(lModeHi-lModeLo);

1995

lIntercept := (abs(lModeHi-lModeLo)-(lModeLo))*lSlope ; //make mode lo = 1;

1996

end else begin

1997

DetectMeanStDev (lImg, lVolVox, lMean,lStDev);

1998

if lStDev <>0 then

1999

lSlope := 1/lStDev

2000

else begin

2001

Msg('Warning: StDev = 0!!!!');

2002

lSlope := 1;

2003

end;

2004

lIntercept := (-lMean*lSlope)+2; //mean voxel has intensity of zero

2005

2006

Msg('method for intensity normalization: Mean = 2, StDev = 1');

2007

Msg('raw_Mean = '+floattostr(lMean)+' '+' raw_StDev = '+floattostr(lStDev));

2008

2009

end;

2010

Msg('Slope/Intercept ' +floattostr(lSlope)+'..'+floattostr(lIntercept));

2011

lHdr.NIFTIhdr.scl_slope := lSlope;

2012

lHdr.NIFTIhdr.scl_inter := lIntercept;

2013

//CopyFileEX(lHdr.HdrFilename,changefileext( lHdr.HdrFilename,'.hdx'));

2014

RenameFile(lHdr.HdrFilename,changefileext( lHdr.HdrFilename,'.hdx'));

2015

//optional - save input

2016

lOutNameMod := ChangeFilePrefixExt(lImageName,'i','.hdr');

2017

NIFTIhdr_SaveHdrImg(lOutNameMod,lHdr.NIFTIhdr,true,not IsNifTiMagic(lHdr.NIFTIhdr),true,lImg,1);

2018

//end optional

2019

NIFTIhdr_SaveHdr(lHdr.HdrFilename,lHdr.NIFTIhdr,true,not IsNifTiMagic(lHdr.NIFTIhdr));

2020

2021

freemem(lImg);

2022

end;

2023

2024

procedure TMainForm.Balance1Click(Sender: TObject);

2025

var

2026

lFilename,lMaskName: string;

2027

lPos: Integer;

2028

begin

2029

MsgClear;

2030

Msg(GetKVers);

2031

if not OpenDialogExecute('Select images for intensity normalization',true,false,kImgFilter) then begin

2032

showmessage('NPM aborted: file selection failed.');

2033

exit;

2034

end; //if not selected

2035

if OpenHdrDlg.Files.Count < 1 then

2036

exit;

2037

lMaskName := '';

2038

for lPos := 1 to OpenHdrDlg.Files.Count do begin

2039

lFilename := OpenHdrDlg.Files[lPos-1];

2040

balance(lFilename,lMaskname,(Sender as TMenuItem).tag);

2041

end;

2042

end;

2043

2044

procedure TMainForm.niiniigz1Click(Sender: TObject);

2045

var

2046

lFilename,lOutname,lPath,lName,lExt: string;

2047

lPos: Integer;

2048

begin

2049

MsgClear;

2050

Msg(GetKVers);

2051

if not OpenDialogExecute('Select images',true,false,kNIIFilter) then begin

2052

showmessage('NPM aborted: file selection failed.');

2053

exit;

2054

end; //if not selected

2055

if OpenHdrDlg.Files.Count < 1 then

2056

exit;

2057

2058

for lPos := 1 to OpenHdrDlg.Files.Count do begin

2059

lFilename := OpenHdrDlg.Files[lPos-1];

2060

FilenameParts(lFilename,lPath,lName,lExt);

2061

lOutname := lPath+lName+'.nii.gz';

2062

msg('Compressing '+ lFilename+' -> '+lOutname);

2063

GZipFile(lFilename, lOutname,false);

2064

end;

2065

msg('Compression completed');

2066

end;

2067

2068

procedure TMainForm.Variance1Click(Sender: TObject);

2069

label

2070

666;

2071

var

2072

lMaskVoxels: integer;

2073

lG: TStrings;

2074

lMaskname: string;

2075

lMaskHdr: TMRIcroHdr;

2076

begin

2077

MsgClear;

2078

Msg(GetKVers);

2079

if not OpenDialogExecute('Select 2 images)',true,true,kImgFilter) then begin

2080

showmessage('NPM aborted: file selection failed.');

2081

exit;

2082

end; //if not selected

2083

lG:= TStringList.Create;

2084

lG.addstrings(OpenHdrDlg.Files);

2085

if lG.count <> 2 then begin

2086

showmessage('You must select exactly two image.');

2087

goto 666;

2088

end;

2089

lMaskname := lG[0];

2090

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

2091

showmessage('Error reading mask.');

2092

goto 666;

2093

end;

2094

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

2095

if not CheckVoxelsGroup(lG,lMaskVoxels) then begin

2096

showmessage('File dimensions differ from mask.');

2097

goto 666;

2098

end;

2099

Msg('Voxels = '+inttostr(lMaskVoxels));

2100

MakeMean(lG,lMaskHdr, odd((Sender as TMenuItem).tag),true);

2101

666:

2102

lG.Free;

2103

end;

2104

2105

procedure BMX;

2106

const

2107

kN = 53;

2108

knNoLesion = 48;

2109

kSymptomRA: array[1..kn] of single =

2110

(4,4.5,2.5,5,4,3.25,0.75,4.5,4.5,0.5,1.625,0,3.5,3,4,2,4.5,5,1.5,5,2.5,5,4,0,2,

2111

1.5,1.75,2.5,5,0,3.25,4.375,0,3.75,0.25,0,2,5,0,0.5,0,2.25,0,2.25,2,0,0.25,0,0,0,0,0,0);

2112

var

2113

lObs: doublep0;

2114

lI: integer;

2115

lBMz,lBMzs,lDF: double;

2116

2117

begin

2118

getmem(lObs,kN * sizeof(double));

2119

for lI := 1 to kN do

2120

lObs^[lI-1] := kSymptomRA[lI];

2121

2122

tBM (kN, knNoLesion, lObs,lBMz,lDF);

2123

//simulate

2124

MainForm.NPMmsg('Generating BM permutation thresholds');

2125

MainForm.Refresh;

2126

genBMsim (kN, lObs);

2127

2128

lBMzs := BMzVal (kN, knNoLesion,lBMz,lDF);

2129

//end simulate

2130

MainForm.NPMmsg('BMsim= '+floattostr(lBMzs)+' '+'BM= '+floattostr(lBMz)+' '+floattostr(lDF) );

2131

freemem(lObs);

2132

end;

2133

2134

(*procedure SX;

2135

var

2136

lVALFilename {,lTemplateName}: string;

2137

lCritPct,lnSubj, lnFactors: integer;

2138

var lSymptomRA: singleP;

2139

var lImageNames: TStrings;

2140

var lCrit: integer; {lBinomial : boolean;}

2141

var lPredictorList: TStringList;

2142

2143

begin

2144

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

2145

lVALFilename := 'c:\RT_norm.val';//MainForm.OpenHdrDlg.Filename;

2146

lPredictorList := TStringList.Create;

2147

GetValCore ( lVALFilename, lnSubj, lnFactors, lSymptomRA, lImageNames, lCrit,lCritPct{,binom},lPredictorList);

2148

lImageNames.Free;

2149

lPredictorList.Free;

2150

end;*)

2151

2152

(*procedure ComputeR;

2153

var

2154

lStr: string;

2155

lT,lDF: double;

2156

begin

2157

inputquery('Enter value','Enter T score',lStr);

2158

lT := strtofloat(lStr);

2159

inputquery('Enter value','Enter DF',lStr);

2160

lDF := strtofloat(lStr);

2161

showmessage('The coresponding correlation Z score for t('+floattostr(lDF)+')='+floattostr(lT) +' is '+floattostr(TtoZ(lT,lDF) ) );

2162

//showmessage('The coresponding correlation R score for t('+floattostr(lDF)+')='+floattostr(lT) +' is '+floattostr(TtoR(lT,lDF) ) );

2163

end;

2164

2165

function Log10x (lLogP: double): double;

2166

begin

2167

result := -Log10(lLogP);

2168

fx(result);

2169

end;

2170

2171

2172

procedure LogPtoZ (lLogP: double);

2173

var

2174

lD,lZ: double;

2175

begin

2176

///lD := Log10(lLogp);

2177

lD := Power(10,-lLogP);

2178

lZ := pNormalInv(lD);

2179

fx(lD,lZ);

2180

end; *)

2181

2182

procedure TMainForm.About1Click(Sender: TObject);

2183

begin

2184

//Masked1Click(nil); exit;

2185

//LogPtoZ (Log10x(0.02));

2186

//LogPtoZ(1.699);

2187

//ComputeR;

2188

showmessage(GetkVers );

2189

end;

2190

2191

procedure TMainForm.Design1Click(Sender: TObject);

2192

begin

2193

{$IFDEF SPREADSHEET} SpreadForm.Show; {$ELSE} Showmessage('Spreadsheet not yet supported on the Operating System');{$ENDIF}

2194

end;

2195

2196

procedure TMainForm.StrToMemo(lStr: String);

2197

var

2198

lLen,lPos: integer;

2199

lOutStr: string;

2200

begin

2201

lLen := length(lStr);

2202

if lLen < 1 then exit;

2203

lOutStr := '';

2204

for lPos := 1 to lLen do begin

2205

if lStr[lPos] = kCR then begin

2206

Msg(lOutStr);

2207

lOutStr := '';

2208

end else

2209

lOutStr := lOutStr + lStr[lPos];

2210

end;

2211

if lOutStr <> '' then

2212

Msg(lOutStr);

2213

end;

2214

2215

2216

procedure TMainForm.PhysiologicalArtifactCorrection1Click(Sender: TObject);

2217

var

2218

lInImgName,lPulsFile,lRespFile,lOutImgName,lStr: string;

2219

l4Ddata: singlep;

2220

lHdr: TMRIcroHdr;

2221

lDim,lImgVox: integer;

2222

lOutHdr: TniftiHdr;

2223

begin

2224

if not OpenDialogExecute('Select file with pulse onsets',false,false,'Siemens physio |*.puls|3-column text |*.txt') then

2225

exit;

2226

lPulsFile := OpenHdrDlg.Filename;

2227

if UpCaseExt(lPulsFile) = '.PULS' then

2228

lRespFile := changefileext(lPulsFile,'.resp')

2229

else begin //text input

2230

if not OpenDialogExecute('Select file with respiration onsets',false,false,'3-column text |*.txt') then

2231

lRespFile := ''

2232

else

2233

lRespFile := OpenHdrDlg.Filename;

2234

end;

2235

if not OpenDialogExecute('Select 4D motion corrected data ',false,false,kImgFilter) then

2236

exit;

2237

lInImgName := OpenHdrDlg.Filename;

2238

if not NIFTIhdr_LoadHdr(lInImgName,lHdr) then begin

2239

showmessage('Error reading image header.');

2240

exit;

2241

end;

2242

for lDim := 1 to 4 do

2243

if lHdr.NIFTIhdr.Dim[lDim] < 4 then begin

2244

Showmessage('You need to select a 4D image with at least 4 voxels/images in each dimension.');

2245

exit;

2246

end;

2247

lImgVox := lHdr.NIFTIhdr.Dim[1]*lHdr.NIFTIhdr.Dim[2]*lHdr.NIFTIhdr.Dim[3];

2248

lDim := lImgVox*lHdr.NIFTIhdr.Dim[4];

2249

MsgClear;

2250

Msg(kVers);

2251

Msg('Physiological Artifact Removal Tool started = ' +TimeToStr(Now));

2252

Msg('Assuming continuous fMRI ascending acquisition with TR = '+realtostr(lHdr.NIFTIhdr.PixDim[4],4)+'sec');

2253

MainForm.refresh;

2254

getmem(l4Ddata,lDim*sizeof(single));

2255

if not LoadImg(lHdr.ImgFileName, l4Ddata, 1, lDim,round(lHdr.NIFTIhdr.vox_offset),1,lHdr.NIFTIhdr.datatype,lImgVox) then begin

2256

Showmessage('Unable to load data');

2257

freemem(l4Ddata);

2258

exit;

2259

end;

2260

lStr := ApplyPART( lPulsFile,l4Ddata,40, lImgVox,lHdr.NIFTIhdr.Dim[3], lHdr.NIFTIhdr.Dim[4], lHdr.NIFTIhdr.PixDim[4]);

2261

if lStr = '' then begin

2262

Showmessage('Unable to apply physio file. Physiological correction is being aborted.');

2263

exit;

2264

end;

2265

StrToMemo (lStr);

2266

if (lRespFile <> '') and (fileexists(lRespFile)) then begin

2267

lStr := ApplyPART( lRespFile,l4Ddata,20, lImgVox,lHdr.NIFTIhdr.Dim[3], lHdr.NIFTIhdr.Dim[4], lHdr.NIFTIhdr.PixDim[4]);

2268

StrToMemo (lStr);

2269

if lStr = '' then begin

2270

Showmessage('Unable to read Respiration file. Hysiological correction is being aborted.');

2271

exit;

2272

end;

2273

2274

end;

2275

2276

MakeHdr (lHdr.NIFTIhdr,lOutHdr);

2277

Msg('Input = ' +lInImgName);

2278

lOutImgName := ChangeFilePrefixExt(lInImgName,'i','.hdr');

2279

NIFTIhdr_SaveHdrImg(lOutImgName,lOutHdr,true,not IsNifTiMagic(lHdr.NIFTIhdr),true,l4Ddata,lHdr.NIFTIhdr.Dim[4]);

2280

Msg('Output = ' +lOutImgName);

2281

Msg('Physiological Artifact Removal Tool finished = ' +TimeToStr(Now));

2282

lOutImgName := ChangeFilePostfixExt(lOutImgName,'Notes','.txt');

2283

MsgSave(lOutImgName);

2284

freemem(l4Ddata);

2285

end;

2286

2287

function ChangeName (lInName: string): string;

2288

var

2289

lPath,lName,lExt: string;

2290

begin

2291

//lInName:= 'c:\vbm\ds\123';

2292

FilenameParts (lInName, lPath,lName,lExt);

2293

//showmessage(lPath+'*'+lName+'*'+lExt);

2294

if length(lName) > 0 then

2295

lName[1] := 'e'

2296

else

2297

lName := 'Unable to convert '+lInName;

2298

result := lPath+lName+lExt;

2299

end;

2300

2301

function Add2ndScans(var lImageNames: TStrings): boolean;

2302

var

2303

lnSubj,lSubj: integer;

2304

lFilename: string;

2305

begin

2306

result := false;

2307

lnSubj :=lImageNames.Count;

2308

if lnSubj < 1 then

2309

exit;

2310

for lSubj := 1 to lnSubj do begin

2311

lFilename := ChangeName(lImageNames[lSubj-1]);

2312

if not (fileexists4D(lFilename)) then begin

2313

showmessage('Unable to find a file named '+ lFilename);

2314

exit;

2315

end;

2316

lImageNames.add(lFilename);

2317

end;

2318

result := true;

2319

end;

2320

2321

function ReadPairedFilenames(var lImageNames: TStrings): boolean;

2322

var

2323

lLen,lPos: integer;

2324

lFilenames,lF1,lF2: string;

2325

lImageNames2: TStrings;

2326

lF: TextFile;

2327

begin

2328

result := false;

2329

Showmessage('Please select a text file with the image names. '+kCR+

2330

'Each line of the file should specify the control and experimental filenames, separated by an *'+kCR+

2331

'C:\vbmdata\c1.nii.gz*C:\vbmdata\e1.nii.gz'+kCR +

2332

'C:\vbmdata\c2.nii.gz*C:\vbmdata\e2.nii.gz'+kCR+

2333

'C:\vbmdata\c3.nii.gz*C:\vbmdata\e3.nii.gz'+kCR+

2334

'...' );

2335

if not MainForm.OpenDialogExecute('Select asterix separated filenames ',false,false,kTxtFilter) then

2336

exit;

2337

lImageNames2:= TStringList.Create; //not sure why TStrings.Create does not work???

2338

//xxx

2339

assignfile(lF,MainForm.OpenHdrDlg.FileName );

2340

FileMode := 0; //read only

2341

reset(lF);

2342

while not EOF(lF) do begin

2343

readln(lF,lFilenames);

2344

lLen := length(lFilenames);

2345

2346

if lLen > 0 then begin

2347

lF1:= '';

2348

lF2 := '';

2349

lPos := 1;

2350

while (lPos <= lLen) and (lFilenames[lPos] <> '*') do begin

2351

lF1 := lF1 + lFilenames[lPos];

2352

inc(lPos);

2353

end;

2354

inc(lPos);

2355

while (lPos <= lLen) do begin

2356

lF2 := lF2 + lFilenames[lPos];

2357

inc(lPos);

2358

end;

2359

if (length(lF1) > 0) and (length(lF2)>0) then begin

2360

if Fileexists4D(lF1) then begin

2361

if Fileexists4D(lF2) then begin

2362

lImageNames.add(lF1);

2363

lImageNames2.add(lF2);

2364

end else //F2exists

2365

showmessage('Can not find image '+lF2);

2366

end else //F1 exists

2367

showmessage('Can not find image '+lF1);

2368

end;

2369

end;//len>0

2370

end; //while not EOF

2371

closefile(lF);

2372

FileMode := 2; //read/write

2373

if (lImageNames.count > 0) and (lImageNames2.count = lImageNames.count) then begin

2374

lImageNames.AddStrings(lImageNames2);

2375

2376

result := true;

2377

end;

2378

lImageNames2.Free;

2379

end;

2380

2381

function AddNumStr(var X : PMatrix; var lNumStr: string; lRow,lCol: integer):boolean;

2382

var

2383

lTempFloat: double;

2384

begin

2385

2386

result := false;

2387

if (lNumStr = '') or (lRow < 1) or (lCol < 1) then exit;

2388

try

2389

lTempFloat := strtofloat(lNumStr);

2390

except

2391

on EConvertError do begin

2392

showmessage('Empty cells? Error reading TXT file row:'+inttostr(lRow)+' col:'+inttostr(lCol)+' - Unable to convert the string '+lNumStr+' to a number');

2393

exit;

2394

end;

2395

end;

2396

//fx(lRow,lCol,lTempFloat);

2397

X^[lCol]^[lRow] := lTempFloat;

2398

lNumStr := '';

2399

result := true;

2400

end;

2401

2402

{$DEFINE notRTEST}

2403

function ReadPairedFilenamesReg(var lImageNames: TStrings; var X : PMatrix; var lnAdditionalFactors: integer): boolean;

2404

var

2405

lLen,lPos,lSep,lMaxSep,lLine: integer;

2406

lFilenames,lF1,lF2,lNumStr: string;

2407

lImageNames2: TStrings;

2408

lF: TextFile;

2409

begin

2410

result := false;

2411

{$IFDEF RTEST}

2412

MainForm.OpenHdrDlg.FileName := 'c:\twins\dataplus.txt';

2413

{$ELSE}

2414

Showmessage('Please select a text file with the image names. '+kCR+

2415

'Each line of the file should specify the control and experimental filenames, separated by an *'+kCR+

2416

'C:\vbmdata\c1.nii.gz*C:\vbmdata\e1.nii.gz'+kCR +

2417

'C:\vbmdata\c2.nii.gz*C:\vbmdata\e2.nii.gz'+kCR+

2418

'C:\vbmdata\c3.nii.gz*C:\vbmdata\e3.nii.gz'+kCR+

2419

'...' );

2420

if not MainForm.OpenDialogExecute('Select asterix separated filenames ',false,false,kTxtFilter) then

2421

exit;

2422

{$ENDIF}

2423

2424

lImageNames2:= TStringList.Create; //not sure why TStrings.Create does not work???

2425

//xxx

2426

assignfile(lF,MainForm.OpenHdrDlg.FileName );

2427

FileMode := 0; //read only

2428

reset(lF);

2429

while not EOF(lF) do begin

2430

readln(lF,lFilenames);

2431

lLen := length(lFilenames);

2432

2433

if lLen > 0 then begin

2434

lF1:= '';

2435

lF2 := '';

2436

lPos := 1;

2437

while (lPos <= lLen) and (lFilenames[lPos] <> '*') do begin

2438

lF1 := lF1 + lFilenames[lPos];

2439

inc(lPos);

2440

end;

2441

inc(lPos);

2442

while (lPos <= lLen) and (lFilenames[lPos] <> '*') do begin

2443

lF2 := lF2 + lFilenames[lPos];

2444

inc(lPos);

2445

end;

2446

if (length(lF1) > 0) and (length(lF2)>0) then begin

2447

if Fileexists4D(lF1) then begin

2448

if Fileexists4D(lF2) then begin

2449

lImageNames.add(lF1);

2450

lImageNames2.add(lF2);

2451

end else //F2exists

2452

showmessage('Can not find image '+lF2);

2453

end else //F1 exists

2454

showmessage('Can not find image '+lF1);

2455

end;

2456

end;//len>0

2457

end; //while not EOF

2458

2459

//fx(lImageNames.count);

2460

//next - count additional factors

2461

lnAdditionalFactors := 0;

2462

reset(lF);

2463

lMaxSep := 0;

2464

while not EOF(lF) do begin

2465

readln(lF,lFilenames);

2466

lLen := length(lFilenames);

2467

lSep := 0;

2468

if lLen > 0 then begin

2469

for lPos := 1 to lLen do

2470

if lFilenames[lPos] = '*' then

2471

inc(lSep)

2472

end;//len>0

2473

if lSep > lMaxSep then

2474

lMaxSep := lSep;

2475

end; //while not EOF

2476

if (lMaxSep > 1) and (lImageNames2.count > 1) then begin //additional factors present

2477

//final pas - load additional factors

2478

lnAdditionalFactors := lMaxSep - 1;

2479

2480

DimMatrix(X, lnAdditionalFactors, lImageNames2.count);

2481

reset(lF);

2482

lLine := 0;

2483

while not EOF(lF) do begin

2484

readln(lF,lFilenames);

2485

lLen := length(lFilenames);

2486

lSep := 0;

2487

2488

if lLen > 0 then begin

2489

inc(lLine);

2490

lPos := 1;

2491

lNumStr := '';

2492

while lPos <= lLen do begin

2493

if (lFilenames[lPos] = '*') then begin

2494

AddNumStr(X,lNumStr,lLine,lSep-1);

2495

inc(lSep);

2496

end else if (lSep >= 2) and (not (lFilenames[lPos] in [#10,#13,#9]) ) then begin

2497

lNumStr := lNumStr+lFilenames[lPos];

2498

//showmessage(lNumStr);

2499

end;

2500

inc(lPos);

2501

end; //while not EOLN

2502

AddNumStr(X,lNumStr,lLine,lSep-1);

2503

end;//len>0

2504

end; //while not EOF

2505

//next - read final line of unterminated string...

2506

end;//maxsepa > 1

2507

2508

2509

//2nd pass vals

2510

closefile(lF);

2511

FileMode := 2; //read/write

2512

if (lImageNames.count > 0) and (lImageNames2.count = lImageNames.count) then begin

2513

lImageNames.AddStrings(lImageNames2);

2514

2515

result := true;

2516

end;

2517

lImageNames2.Free;

2518

result := true;

2519

end;

2520

2521

procedure TMainForm.DualImageCorrelation1Click(Sender: TObject);

2522

label

2523

666;

2524

var

2525

lnSubj,lSubj,lMaskVoxels,lnAdditionalFactors,lI: integer;

2526

lImageNames: TStrings;

2527

X: PMatrix;

2528

lMaskname,lStr,lOutName: string;

2529

lMaskHdr: TMRIcroHdr;

2530

begin

2531

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

2532

MsgClear;

2533

Msg(kVers);

2534

{$IFDEF RTEST}

2535

OpenHdrDlg.FileName := 'c:\twins\aameanMean.hdr';

2536

{$ELSE}

2537

Msg('Dual-image Linear Regression [Weighted Least Squares]');

2538

if not OpenDialogExecute('Select brain mask ',false,false,kImgFilter) then begin

2539

showmessage('NPM aborted: mask selection failed.');

2540

goto 666;

2541

end; //if not selected

2542

{$ENDIF}

2543

2544

lMaskname := OpenHdrDlg.Filename;

2545

2546

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

2547

showmessage('Error reading Mask image.');

2548

goto 666;

2549

end;

2550

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

2551

if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file} then begin

2552

showmessage('Mask file size too small.');

2553

goto 666;

2554

end;

2555

if not ReadPairedFilenamesReg(lImageNames,X,lnAdditionalFactors) then exit;

2556

lnSubj :=lImageNames.Count div 2;

2557

2558

//fx(lnAdditionalFactors);

2559

//show matrix

2560

//MsgStrings (lImageNames);

2561

Msg ('n Subjects = '+inttostr(lnSubj));

2562

for lSubj := 0 to (lnSubj-1) do begin

2563

lStr := lImageNames[lSubj]+' <-> '+lImageNames[lSubj+lnSubj];

2564

if lnAdditionalFactors > 0 then

2565

for lI := 1 to lnAdditionalFactors do

2566

lStr := lStr+','+floattostr(X^[lI]^[lSubj+1]);

2567

2568

2569

Msg(lStr);

2570

end;

2571

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

2572

showmessage('File dimensions differ from mask.');

2573

goto 666;

2574

end;

2575

2576

2577

Msg('Mask = '+lMaskname);

2578

Msg('Total voxels = '+inttostr(lMaskVoxels));

2579

Msg('Number of observations = '+inttostr(lnSubj));

2580

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

2581

showmessage('File dimensions differ from mask.');

2582

goto 666;

2583

end;

2584

2585

if lnSubj < 5 then begin

2586

showmessage('Paired regression error: Requires at least 5 images per group.');

2587

goto 666;

2588

end;

2589

lOutName := lMaskName;

2590

if not SaveHdrName ('Base Statistical Map', lOutName) then exit;

2591

//showmessage('Unimplemented Regress');//

2592

Regress2NPMAnalyze (lImageNames, lMaskHdr, lOutname,X,lnAdditionalFactors);

2593

if lnAdditionalFactors > 1 then

2594

DelMatrix(X, lnAdditionalFactors, lnSubj);

2595

666:

2596

lImageNames.Free;

2597

end;

2598

2599

procedure TMainForm.LesionBtnClick(Sender: TObject);

2600

label

2601

666;

2602

var

2603

lPrefs: TLDMPrefs ;

2604

begin

2605

lPrefs.NULP := gNULP;

2606

if (1= (Sender as tMenuItem).tag) then begin //continuous

2607

lPrefs.BMtest := BMmenu.checked;

2608

lPrefs.Ttest := ttestmenu.checked;

2609

if (not lPrefs.BMtest) and (not lPrefs.ttest) then

2610

lPrefs.ttest := true;

2611

lPrefs.Ltest:= false;

2612

end else begin //binomial

2613

lPrefs.BMtest := false;

2614

lPrefs.Ttest := false;

2615

lPrefs.Ltest:= true;

2616

end;

2617

lPrefs.CritPct := -1;

2618

lPrefs.nPermute := ReadPermute;

2619

lPrefs.Run := 0;{0 except for montecarlo}

2620

{if (not lPrefs.Ltest) and (not lPrefs.Ttest) and (not lPrefs.BMtest) then begin

2621

Showmessage('Error: you need to compute at least on test [options/test menu]');

2622

exit;

2623

end; code above defaults to t-test}

2624

if not MainForm.OpenDialogExecute('Select MRIcron VAL file',false,false,'MRIcron VAL (*.val)|*.val') then begin

2625

showmessage('NPM aborted: VAL file selection failed.');

2626

exit;

2627

end; //if not selected

2628

lPrefs.VALFilename := MainForm.OpenHdrDlg.Filename;

2629

lPrefs.OutName := ExtractFileDirWithPathDelim(lPrefs.VALFilename)+'results';

2630

lPrefs.OutName := lPrefs.OutName+'.nii.gz';

2631

SaveHdrDlg.Filename := lPrefs.Outname;

2632

if not SaveHdrName ('Base Statistical Map', lPrefs.OutName) then exit;

2633

DoLesion (lPrefs); //Prefs.pas

2634

end;

2635

2636

2637

procedure TMainForm.FormShow(Sender: TObject);

2638

begin

2639

MsgClear;

2640

Msg(GetkVers);

2641

{$IFNDEF UNIX} {GUILaunch;}{$ENDIF}

2642

LongTimeFormat := 'YYYY-MMM-DD hh:nn:ss'; //delphi TimeToStr

2643

ShortTimeFormat := 'YYYY-MMM-DD hh:nn:ss'; //freepascal TimeToStr

2644

//stax;

2645

//MakeGZ;

2646

//{x$IFNDEF UNIX} Threads1.visible := false;{x$ENDIF}//Windows can read actual CPU count

2647

//TestMultReg;

2648

//SingleRegressClick(nil);

2649

//DualImageCorrelation1Click(nil);

2650

//AnaCOM1Click(nil);

2651

//msg(floattostr(pNormalInv(1/20000)));

2652

//msg(floattostr(pNormalInv(0.01667)));

2653

//msg(floattostr(pNormalInv(0.05/51700)));

2654

//msg(floattostr(0.05/pNormal(4.76 )));

2655

{$IFNDEF UNIX} ReadParamStr; {$ENDIF}

2656

//FX(rocA(0.2,0.4),AUC(0.7,0.4),rocA(0.4,0.7),AUC(0.4,0.7) );

2657

//LesionX;

2658

{$IFDEF benchmark}

2659

MonteCarloSimulation1.visible := true;

2660

// LesionMonteCarlo (false,true,true);

2661

{$ENDIF}

2662

end;

2663

2664

2665

procedure TMainForm.PairedTMenuClick(Sender: TObject);

2666

label

2667

666;

2668

var

2669

lnSubj,lSubj,lMaskVoxels: integer;

2670

lImageNames: TStrings;

2671

lMaskname,lStr,lOutName: string;

2672

lMaskHdr: TMRIcroHdr;

2673

begin

2674

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

2675

MsgClear;

2676

Msg(kVers);

2677

Msg('Paired T-test [Repeated Measures]');

2678

if not OpenDialogExecute('Select brain mask ',false,false,kImgFilter) then begin

2679

showmessage('NPM aborted: mask selection failed.');

2680

goto 666;

2681

end; //if not selected

2682

//OpenHdrDlg.FileName := 'c:\vbmdata\mask50.nii.gz';

2683

lMaskname := OpenHdrDlg.Filename;

2684

2685

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

2686

showmessage('Error reading Mask image.');

2687

goto 666;

2688

end;

2689

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

2690

if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file} then begin

2691

showmessage('Mask file size too small.');

2692

goto 666;

2693

end;

2694

if not ReadPairedFilenames(lImageNames) then exit;

2695

lnSubj :=lImageNames.Count div 2;

2696

2697

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

2698

showmessage('File dimensions differ from mask.');

2699

goto 666;

2700

end;

2701

2702

2703

Msg('Mask = '+lMaskname);

2704

Msg('Total voxels = '+inttostr(lMaskVoxels));

2705

Msg('Number of observations = '+inttostr(lnSubj));

2706

Msg('Degrees of Freedom = '+inttostr(lnSubj-1));

2707

2708

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

2709

showmessage('File dimensions differ from mask.');

2710

goto 666;

2711

end;

2712

//show matrix

2713

//MsgStrings (lImageNames);

2714

Msg ('n Subjects = '+inttostr(lnSubj));

2715

lStr := 'Image,';

2716

for lSubj := 0 to (lnSubj-1) do

2717

Msg(lImageNames[lSubj]+' <-> '+lImageNames[lSubj+lnSubj]);

2718

if lnSubj < 4 then begin

2719

showmessage('Paired t-test error: Requires at least 4 images per group.');

2720

goto 666;

2721

end;

2722

lOutName := lMaskName;

2723

//if not SaveHdrName ('Base Statistical Map', lOutName) then exit;

2724

NPMAnalyzePaired (lImageNames, lMaskHdr, lMaskVoxels);

2725

//Regress2NPMAnalyze (lImageNames, lMaskHdr, lOutname);

2726

666:

2727

lImageNames.Free;

2728

end;

2729

2730

function TMainForm.NPMzscore (var lImages: TStrings; var lMnHdr,lStDevHdr: TMRIcroHdr): boolean;

2731

label

2732

667;

2733

var

2734

lOutNameMod: string;

2735

lMnImg,lStDevImg,lSubjImg,lOutImg: SingleP;

2736

lVal: single;

2737

lSubj,lPos,lVolVox: integer;

2738

lStatHdr: TNIfTIhdr;

2739

begin

2740

result := false;

2741

//lOutName := lMnHdr.ImgFileName;

2742

//if not SaveHdrName ('Statistical Map', lOutName) then exit;

2743

Msg('Analysis began = ' +TimeToStr(Now));

2744

lVolVox := lMnHdr.NIFTIhdr.dim[1]*lMnHdr.NIFTIhdr.dim[2]* lMnHdr.NIFTIhdr.dim[3];

2745

if (lVolVox < 1) then goto 667;

2746

//load mask

2747

for lPos := 0 to lImages.Count do

2748

if gScaleRA[lPos] = 0 then

2749

gScaleRA[lPos] := 1;

2750

if gScaleRA[kMaxImages] = 0 then

2751

gScaleRA[kMaxImages] := 1;

2752

2753

getmem(lMnImg,lVolVox*sizeof(single));

2754

if not LoadImg(lMnHdr.ImgFileName, lMnImg, 1, lVolVox,round(gOffsetRA[0]),1,lMnHdr.NIFTIhdr.datatype,lVolVox) then begin

2755

Msg('Unable to load mean ' +lMnHdr.ImgFileName);

2756

goto 667;

2757

end;

2758

//load StDev

2759

getmem(lStDevImg,lVolVox*sizeof(single));

2760

if not LoadImg(lStDevHdr.ImgFileName, lStDevImg, 1, lVolVox,round(gOffsetRA[kMaxImages]),1,lStDevHdr.NIFTIhdr.datatype,lVolVox) then begin

2761

Msg('Unable to load StDev ' +lStDevHdr.ImgFileName);

2762

goto 667;

2763

end;

2764

getmem(lOutImg,lVolVox* sizeof(single));

2765

for lPos := 1 to lVolVox do begin

2766

lMnImg^[lPos] := (gScaleRA[0]*lMnImg^[lPos])+gInterceptRA[0];

2767

lStDevImg^[lPos] := (gScaleRA[kMaxImages]*lStDevImg^[lPos])+gInterceptRA[kMaxImages];

2768

if lStDevImg^[lPos] = 0 then

2769

lOutImg^[lPos] := 0;

2770

end;

2771

getmem(lSubjImg,lVolVox* sizeof(single));

2772

for lSubj := 1 to lImages.Count do begin

2773

ProgressBar1.Position := round((lSubj/lImages.Count)*100);

2774

Msg( lImages.Strings[lSubj-1]);

2775

showmessage(inttostr(round(gOffsetRA[lSubj])));

2776

LoadImg(lImages.Strings[lSubj-1], lSubjImg, 1, lVolVox,round(gOffsetRA[lSubj]),1,gDataTypeRA[lSubj],lVolVox);

2777

for lPos := 1 to lVolVox do begin

2778

if lStDevImg^[lPos] <> 0 then begin

2779

lVal := (gScaleRA[lSubj]*lSubjImg^[lPos])+gInterceptRA[lSubj];

2780

lOutImg^[lPos] := (lVal-lMnImg^[lPos])/lStDevImg^[lPos];

2781

end; //for each voxel with variance

2782

end; //for each voxel

2783

lOutNameMod := ChangeFilePostfixExt(lImages.Strings[lSubj-1],'Z','.hdr');

2784

MakeStatHdr (lMnHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lVolVox,kNIFTI_INTENT_ZSCORE,inttostr(lVolVox) );

2785

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMnHdr.NIFTIhdr),true,lOutImg,1);

2786

end; //for each subj

2787

freemem(lSubjImg);

2788

freemem(lOutImg);

2789

freemem(lMnImg);

2790

freemem(lStDevImg);

2791

Msg('Analysis finished = ' +TimeToStr(Now));

2792

ProgressBar1.Position := 0;

2793

result := true;

2794

exit;

2795

667: //you only get here if you aborted ... free memory and report error

2796

if lVolVox > 1 then freemem(lMnImg);

2797

Msg('Unable to complete analysis.');

2798

ProgressBar1.Position := 0;

2799

end;

2800

2801

2802

procedure TMainForm.SingleSubjectZScores1Click(Sender: TObject);

2803

label

2804

666;

2805

var

2806

lnSubj,lMnVoxels: integer;

2807

lG: TStrings;

2808

lMn,lStDev: string;

2809

lMnHdr,lStDevHdr: TMRIcroHdr;

2810

begin

2811

if (not ttestmenu.checked) and (not BMmenu.checked) then begin

2812

Showmessage('Error: you need to compute at least on test [options/test menu]');

2813

exit;

2814

end;

2815

MsgClear;

2816

Msg(kVers);

2817

Msg('Threads: '+inttostr(gnCPUThreads));

2818

if not OpenDialogExecute('Select mean image ',false,false,kImgFilter) then begin

2819

showmessage('NPM aborted: mean selection failed.');

2820

exit;

2821

end; //if not selected

2822

lMn := OpenHdrDlg.Filename;

2823

if not NIFTIhdr_LoadHdr(lMn,lMnHdr) then begin

2824

showmessage('Error reading mask.');

2825

exit;

2826

end;

2827

lMnVoxels := ComputeImageDataBytes8bpp(lMnHdr);

2828

if (lMnVoxels < 2) or (not CheckVoxels(lMn,lMnVoxels,0)){make sure there is uncompressed .img file} then begin

2829

showmessage('Mean file size too small.');

2830

exit;

2831

end;

2832

2833

if not OpenDialogExecute('Select StDev image ',false,false,kImgFilter) then begin

2834

showmessage('NPM aborted: StDev selection failed.');

2835

exit;

2836

end; //if not selected

2837

lStDev := OpenHdrDlg.Filename;

2838

if not NIFTIhdr_LoadHdr(lStDev,lStDevHdr) then begin

2839

showmessage('Error reading StDev.');

2840

exit;

2841

end;

2842

if not CheckVoxels(lStDev, lMnVoxels,kMaxImages) then begin

2843

showmessage('Error Mean and StDev must have same size.');

2844

exit;

2845

end;

2846

Msg('Mean name = '+ lMn);

2847

Msg('Total voxels = '+inttostr(lMnVoxels));

2848

//next, get 1st group

2849

if not OpenDialogExecute('Select postive group (Z scores positive if this group is brighter)',true,false,kImgFilter) then begin

2850

showmessage('NPM aborted: file selection failed.');

2851

exit;

2852

end; //if not selected

2853

lG:= TStringList.Create; //not sure why TStrings.Create does not work???

2854

lG.addstrings(OpenHdrDlg.Files);

2855

lnSubj :=OpenHdrDlg.Files.Count;

2856

Msg('Subjects= '+inttostr(lnSubj));

2857

if not CheckVoxelsGroup(lG,lMnVoxels) then begin

2858

showmessage('File dimensions differ from mask.');

2859

goto 666;

2860

end;

2861

NPMzscore (lG, lMnHdr,lStDevHdr);

2862

//NPMAnalyze(lG,lMnHdr,lMnVoxels,lnSubj);

2863

666:

2864

lG.Free;

2865

end;

2866

2867

procedure TMainForm.MultipleRegressClick(Sender: TObject);

2868

label

2869

666;

2870

var

2871

lnFactors,lnSubj,lMaskVoxels,lRow,lCol: integer;

2872

lImageNames: TStrings;

2873

lPredictorList: TStringList;

2874

lTemp4D,lMaskname,lStr,lOutName: string;

2875

lMaskHdr: TMRIcroHdr;

2876

X : PMatrix;

2877

begin

2878

{$IFDEF FPC}

2879

showmessage('Regression routines not extensively tested: you may want to use the Windows compilation.');

2880

{$ENDIF}

2881

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

2882

lPredictorList := TStringList.Create;

2883

Memo1.Lines.Clear;

2884

Memo1.Lines.Add(kVers);

2885

Memo1.Lines.Add('Multiple Linear Regression [Weighted Least Squares]');

2886

if not GetValReg(lnSubj,lnFactors,X,lImageNames,lPredictorList) then

2887

goto 666;

2888

lTemp4D := CreateDecompressed4D(lImageNames);

2889

if not OpenDialogExecute('Select brain mask ',false,false,kImgFilter) then begin

2890

showmessage('NPM aborted: mask selection failed.');

2891

goto 666;

2892

end; //if not selected

2893

lMaskname := OpenHdrDlg.Filename;

2894

//lMaskname := lImageNames[0];

2895

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

2896

showmessage('Error reading 1st image.');

2897

goto 666;

2898

end;

2899

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

2900

if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file} then begin

2901

showmessage('Mask file size too small.');

2902

goto 666;

2903

end;

2904

Memo1.Lines.Add('Mask = '+lMaskname);

2905

Memo1.Lines.Add('Total voxels = '+inttostr(lMaskVoxels));

2906

Memo1.Lines.Add('Number of observations = '+inttostr(lnSubj));

2907

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

2908

showmessage('File dimensions differ from mask.');

2909

goto 666;

2910

end;

2911

//show matrix

2912

lStr := 'Image,';

2913

for lCol := 1 to lnFactors do

2914

lStr := lStr + lPredictorList.Strings[lCol-1]+', ';

2915

MainForm.Memo1.Lines.Add(lStr);

2916

for lRow := 1 to lnSubj do begin

2917

lStr := lImageNames[lRow-1]+',';

2918

for lCol := 1 to lnFactors do

2919

lStr := lStr + floattostr(X^[lCol]^[lRow])+', ';

2920

MainForm.Memo1.Lines.Add(lStr);

2921

end;

2922

lOutName := lMaskName;

2923

if not SaveHdrName ('Base Statistical Map', lOutName) then exit;

2924

ARegressNPMAnalyze(lImageNames,lMaskHdr,X,lnFactors,lPredictorList,lOutName);

2925

DelMatrix(X, lnFactors, lnSubj);

2926

666:

2927

lImageNames.Free;

2928

lPredictorList.Free;

2929

DeleteDecompressed4D(lTemp4D);

2930

end;

2931

2932

{$DEFINE notRegTest}

2933

procedure TMainForm.SingleRegressClick(Sender: TObject);

2934

label

2935

666;

2936

var

2937

lnSubj1,lnFactors,lnSubj,lMaskVoxels,lRow,lCol: integer;

2938

2939

lImageNames,lImageNames1: TStrings;

2940

lPredictorList,lPredictorList1: TStringList;

2941

lTemp4D,lMaskname,lOutName: string;

2942

lMaskHdr: TMRIcroHdr;

2943

X,X1 : PMatrix;

2944

begin

2945

{$IFDEF FPC}

2946

showmessage('Regression routines not extensively tested: you may want to use the Windows compilation.');

2947

{$ENDIF}

2948

lTemp4D := '';

2949

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

2950

lPredictorList := TStringList.Create;

2951

lPredictorList1 := TStringList.Create;

2952

if not GetValReg(lnSubj,lnFactors,X,lImageNames,lPredictorList) then

2953

goto 666;

2954

{$IFDEF regtest}

2955

lMaskname := 'C:\Documents and Settings\Chris Rorden\Desktop\npmdata\npmdata\amask50.nii.gz';

2956

{$ELSE}

2957

if not OpenDialogExecute('Select brain mask ',false,false,kImgFilter) then begin

2958

showmessage('NPM aborted: mask selection failed.');

2959

goto 666;

2960

end; //if not selected

2961

lMaskname := OpenHdrDlg.Filename;

2962

{$ENDIF}

2963

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

2964

showmessage('Error reading 1st image.');

2965

goto 666;

2966

end;

2967

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

2968

if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file} then begin

2969

showmessage('Mask file size too small.');

2970

goto 666;

2971

end;

2972

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

2973

showmessage('File dimensions differ from mask.');

2974

goto 666;

2975

end;

2976

lOutName := lMaskName;

2977

{$IFNDEF regtest}

2978

2979

if not SaveHdrName ('Base Statistical Map', lOutName) then goto 666;

2980

{$ENDIF}

2981

lTemp4D := CreateDecompressed4D(lImageNames);

2982

2983

2984

lImageNames1:= TStringList.Create;

2985

for lCol := 1 to lnFactors do begin

2986

lPredictorList1.Clear;

2987

lPredictorList1.Add(lPredictorList[lCol-1]);

2988

lImageNames1.clear;

2989

for lRow := 1 to lnSubj do

2990

if X^[lCol]^[lRow] <> kNaN then

2991

lImageNames1.Add(lImageNames[lRow-1]);

2992

DimMatrix(X1, 1, lImageNames1.Count);

2993

lnSubj1 := 0;

2994

2995

for lRow := 1 to lnSubj do

2996

if X^[lCol]^[lRow] <> kNaN then begin

2997

inc(lnSubj1);

2998

X1^[1]^[lnSubj1] := X^[lCol]^[lRow];

2999

end;

3000

if lnSubj1 <> lImageNames1.Count then //should be impossible

3001

showmessage('serious error');

3002

Memo1.Lines.Clear;

3003

Memo1.Lines.Add(kVers);

3004

Memo1.Lines.Add('Single Linear Regression [Weighted Least Squares]');

3005

Memo1.Lines.Add('Mask = '+lMaskname);

3006

Memo1.Lines.Add('Total voxels = '+inttostr(lMaskVoxels));

3007

Memo1.Lines.Add('Number of observations = '+inttostr(lnSubj1));

3008

3009

MainForm.Memo1.Lines.Add('Image,'+ lPredictorList1.Strings[0]);

3010

for lRow := 1 to lnSubj1 do

3011

MainForm.Memo1.Lines.Add(lImageNames1[lRow-1]+','+floattostr(X1^[1]^[lRow]) ) ;

3012

ARegressNPMAnalyze(lImageNames1,lMaskHdr,X1,1,lPredictorList1,lOutName);

3013

DelMatrix(X1, 1, lnSubj1);

3014

end;

3015

lImageNames1.Free;

3016

DelMatrix(X, lnFactors, lnSubj);

3017

666:

3018

DeleteDecompressed4D(lTemp4D);

3019

lImageNames.Free;

3020

lPredictorList.Free;

3021

lPredictorList1.Free;

3022

end;

3023

3024

procedure TMainForm.AssociatevalfileswithNPM1Click(Sender: TObject);

3025

begin

3026

{$IFDEF FPC}

3027

//unsupported by FreePascal

3028

{$ELSE}

3029

case MessageDlg('NPM installation:'+kCR+'Do you want .val fiels to automatically open NPM when you double click on their icons?', mtConfirmation,

3030

[mbYes, mbNo], 0) of { produce the message dialog box }

3031

id_No: exit;

3032

end;

3033

registerfiletype(kVALNativeExt,'NPM'{key},'NPM',Application.ExeName+',1');

3034

{$ENDIF}

3035

3036

end;

3037

3038

procedure TMainForm.threadChange(Sender: TObject);

3039

begin

3040

(sender as tmenuitem).checked := true;

3041

ReadThread;

3042

end;

3043

3044

procedure TMainForm.Countlesionoverlaps1Click(Sender: TObject);

3045

label

3046

666;

3047

var

3048

lReps,lMax,lInc,lMaskVoxels,lDefault,lTotal,lPct: integer;

3049

lG: TStrings;

3050

lMaskname: string;

3051

lMaskHdr: TMRIcroHdr;

3052

begin

3053

MsgClear;

3054

Msg(kVers);

3055

if not OpenDialogExecute('Select images to overlap',true,false,kImgFilter) then begin

3056

showmessage('NPM aborted: file selection failed.');

3057

exit;

3058

end; //if not selected

3059

if MainForm.OpenHdrDlg.Files.Count < 2 then begin

3060

lTotal := NIFTIhdr_HdrVolumes(MainForm.OpenHdrDlg.Filename);

3061

if lTotal < 2 then begin

3062

Showmessage('Error: This function is designed to overlay MULTIPLE volumes. You selected less than two images.');

3063

exit;

3064

end;

3065

lG:= TStringList.Create;

3066

for lReps := 1 to lTotal do

3067

lG.Add(MainForm.OpenHdrDlg.Filename+':'+inttostr(lReps) );

3068

end else begin

3069

lG:= TStringList.Create;

3070

lG.addstrings(OpenHdrDlg.Files);

3071

end;

3072

lMaskname := lG[0];

3073

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

3074

showmessage('Error reading mask.');

3075

goto 666;

3076

end;

3077

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

3078

if not CheckVoxelsGroup(lG,lMaskVoxels) then begin

3079

showmessage('File dimensions differ from mask.');

3080

goto 666;

3081

end;

3082

lTotal := lG.Count;

3083

if lTotal > kMaxObs then

3084

lTotal := kMaxObs; //this implemmentation uses 126 bits per voxel - we can not test more than this!

3085

if lTotal > 100 then

3086

lDefault := 100

3087

else

3088

lDefault := lTotal;

3089

lMax := ReadIntForm.GetInt('Enter maximum number of overlaps to test ', 3,lDefault,lTotal);

3090

lDefault := lMax div 10;

3091

if lDefault < 1 then

3092

lDefault := 1;

3093

lInc := ReadIntForm.GetInt('Enter overlap increment (e.g. if 5; then 5, 10, 15...) ', 1,lDefault,lMax);

3094

lReps := ReadIntForm.GetInt('Enter number of times each increment is tested ', 1,10,100);

3095

lPct := ReadIntForm.GetInt('Only include voxels damaged in N% of patients ', 0,5,100);

3096

3097

Msg('Voxels = '+inttostr(lMaskVoxels));

3098

Msg('Scans to permute = '+inttostr(lG.count));

3099

EvaluatePower (lG,lInc,lMax,lReps,lPct);

3100

3101

//MakeMean(lG,lMaskHdr, odd((Sender as TMenuItem).tag),false);

3102

666:

3103

lG.Free;

3104

end;

3105

3106

{$DEFINE SINGLETHREAD}

3107

{$DEFINE NOTHREAD}

3108

3109

function TMainForm.FirthNPMAnalyze (var lImages: TStrings; var lPredictorList: TStringList; var lMaskHdr: TMRIcroHdr; lnCond,lnCrit: integer; var lSymptomRA: SingleP; var lOutName: string): boolean;

3110

label

3111

123,667;

3112

var

3113

lOutNameMod: string;

3114

lPlankImg: bytep;

3115

lOutImgSum : singleP;

3116

lOutImg: SingleRAp;

3117

{$IFDEF SINGLETHREAD}lnCPUThreads,{$ENDIF}

3118

lCond,lPos,lPlank,lThread,lnDeficit: integer;

3119

lTotalMemory,lVolVox,lMinMask,lMaxMask,lnPlanks,lVoxPerPlank,

3120

lThreadStart,lThreadInc,lThreadEnd, lnLesion,lnPermute,

3121

lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lPosPct: int64;

3122

lT, lSum: double;

3123

lObsp: pointer;

3124

lObs: Doublep0;

3125

lStatHdr: TNIfTIhdr;

3126

lFdata: file;

3127

lRanOrderp: pointer;

3128

lRanOrder: Doublep0;

3129

begin

3130

if lnCond < 1 then

3131

exit;

3132

lnPermute := ReadPermute;

3133

if lnPermute > 1 then begin

3134

Msg('NPM does not (yet) support permutation thresholding with Logisitic Regression.');

3135

lnPermute := 0;

3136

end;

3137

{$IFDEF SINGLETHREAD}

3138

lnCPUThreads := gnCPUThreads;

3139

if gnCPUThreads > 1 then

3140

Msg('July 2007 logistic regression will only use 1 thread. You may want to check for a software update');

3141

gnCPUThreads := 1;

3142

{$ENDIF}

3143

Msg('Permutations = ' +IntToStr(lnPermute));

3144

Msg('Logisitic Regression began = ' +TimeToStr(Now));

3145

lTotalMemory := 0;

3146

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

3147

if (lVolVox < 1) then goto 667;

3148

lMinMask := 1;

3149

lMaxMask := lVolVox;

3150

lVoxPerPlank := kPlankSz div lImages.Count div sizeof(byte) ;

3151

if (lVoxPerPlank = 0) then goto 667; //no data

3152

lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;

3153

if (lTotalMemory = 0) then goto 667; //no data

3154

lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;

3155

Msg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));

3156

Msg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));

3157

if (lnPlanks = 1) then

3158

getmem(lPlankImg,lTotalMemory)

3159

else

3160

getmem(lPlankImg,kPlankSz);

3161

lStartVox := lMinMask;

3162

lEndVox := lMinMask-1;

3163

for lPos := 1 to lImages.Count do

3164

if gScaleRA[lPos] = 0 then

3165

gScaleRA[lPos] := 1;

3166

createArray64(lObsp,lObs,lImages.Count);

3167

getmem(lOutImgSum,lVolVox* sizeof(single));

3168

//getmem(lOutImgL,lVolVox* sizeof(single));

3169

getmem(lOutImg,lnCond*sizeof(Singlep));

3170

for lCond := 1 to lnCond do begin

3171

getmem(lOutImg^[lCond],lVolVox* sizeof(single));

3172

for lPos := 1 to lVolVox do

3173

lOutImg^[lCond]^[lPos] := 0;

3174

end;

3175

//InitPermute (lImages.Count, lnPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);

3176

for lPos := 1 to lVolVox do

3177

lOutImgSum^[lPos] := 0;

3178

ClearThreadData(gnCPUThreads,lnPermute) ;

3179

for lPlank := 1 to lnPlanks do begin

3180

ProgressBar1.Position := 1;

3181

Msg('Computing plank = ' +Inttostr(lPlank));

3182

Refresh;

3183

Application.processmessages;

3184

lEndVox := lEndVox + lVoxPerPlank;

3185

if lEndVox > lMaxMask then begin

3186

lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);

3187

lEndVox := lMaxMask;

3188

end;

3189

lPlankImgPos := 1;

3190

for lPos := 1 to lImages.Count do begin

3191

if not LoadImg8(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then

3192

goto 667;

3193

lPlankImgPos := lPlankImgPos + lVoxPerPlank;

3194

end;//for each image

3195

//threading start

3196

lThreadStart := 1;

3197

lThreadInc := lVoxPerPlank div gnCPUThreads;

3198

lThreadEnd := lThreadInc;

3199

Application.processmessages;

3200

{$IFDEF NOTHREAD}

3201

FirthAnalyzeNoThread (lnCond, lnCrit,lnPermute,1,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,lPlankImg,lOutImgSum,lSymptomRA,lOutImg);

3202

//FirthAnalyzeNoThread (lnCond,lnCrit, lnPermute,1,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,lPlankImg,lOutImgSum,lSymptomRA,lOutImg);

3203

{$ELSE}

3204

for lThread := 1 to gnCPUThreads do begin

3205

if lThread = gnCPUThreads then

3206

lThreadEnd := lVoxPerPlank; //avoid integer rounding error

3207

with TFirthThreadStat.Create (ProgressBar1,lnCond,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,lPlankImg,lOutImgSum,lSymptomRA,lOutImg) do

3208

{$IFDEF FPC} OnTerminate := @ThreadDone; {$ELSE}OnTerminate := ThreadDone;{$ENDIF}

3209

inc(gThreadsRunning);

3210

Msg('Thread ' +Inttostr(gThreadsRunning)+' = '+inttostr(lThreadStart)+'..'+inttostr(lThreadEnd));

3211

lThreadStart := lThreadEnd + 1;

3212

lThreadEnd :=lThreadEnd + lThreadInc;

3213

end; //for each thread

3214

3215

repeat

3216

Application.processmessages;

3217

until gThreadsRunning = 0;

3218

{$ENDIF} //THREADED

3219

Application.processmessages;

3220

//showmessage('Threads done');

3221

//threading end

3222

lStartVox := lEndVox + 1;

3223

end;

3224

lnVoxTested := SumThreadDataLite(gnCPUThreads); //not yet lnVoxTested := SumThreadData(gnCPUThreads,lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);

3225

if lnVoxTested < 1 then begin

3226

Msg('**Error: no voxels tested: no regions lesioned in at least '+inttostr(lnCrit)+' patients**');

3227

goto 123;

3228

end;

3229

//next report findings

3230

Msg('Voxels tested = ' +Inttostr(lnVoxTested));

3231

Msg('Only tested voxels with more than '+inttostr(lnCrit)+' lesions');

3232

//Next: save results from permutation thresholding....

3233

reportBonferroni('Std',lnVoxTested);

3234

//next: save data

3235

(*savedata *)

3236

MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);

3237

//save sum map

3238

lOutNameMod := ChangeFilePostfixExt(lOutName,'Sum','.hdr');

3239

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);

3240

3241

MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );

3242

for lCond := 1 to lnCond do begin

3243

reportFDR (lPredictorList[lCond-1]+inttostr(lCond), lVolVox, lnVoxTested, lOutImg^[lCond]);

3244

//reportPermute('L',lnPermute,lPermuteMaxBM, lPermuteMinBM);

3245

lOutNameMod := ChangeFilePostfixExt(lOutName,lPredictorList[lCond-1]+inttostr(lCond),'.hdr');

3246

NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImg^[lCond],1);

3247

end;

3248

123:

3249

//next: free dynamic memory

3250

//FreePermute (lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp);

3251

for lCond := 1 to lnCond do

3252

freemem(lOutImg^[lCond]);

3253

freemem(lOutImg);

3254

3255

freemem(lOutImgSum);

3256

freemem(lObsp);

3257

freemem(lPlankImg);

3258

Msg('Analysis finished = ' +TimeToStr(Now));

3259

lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes','.txt');

3260

MsgSave(lOutNameMod);

3261

3262

ProgressBar1.Position := 0;

3263

{$IFDEF SINGLETHREAD}

3264

gnCPUThreads := lnCPUThreads;

3265

{$ENDIF}

3266

exit;

3267

667: //you only get here if you aborted ... free memory and report error

3268

if lTotalMemory > 1 then freemem(lPlankImg);

3269

Msg('Unable to complete analysis.');

3270

ProgressBar1.Position := 0;

3271

{$IFDEF SINGLETHREAD}

3272

gnCPUThreads := lnCPUThreads;

3273

{$ENDIF}

3274

end;

3275

3276

3277

function ComputeLesionVolume (lImgName: string): integer;

3278

var

3279

lHdr: TMRIcroHdr;

3280

lImg: byteP;

3281

lVolVox,lVox:integer;

3282

begin

3283

result := -1; //error

3284

NIFTIhdr_LoadHdr(lImgName,lHdr);

3285

lVolVox := lHdr.NIFTIhdr.dim[1]*lHdr.NIFTIhdr.dim[2]* lHdr.NIFTIhdr.dim[3];

3286

getmem(lImg,lVolVox*sizeof(byte));

3287

if not LoadImg8(lImgName, lImg, 1, lVolVox,round(lHdr.NIFTIhdr.vox_offset),1,lHdr.NIFTIhdr.datatype,lVolVox) then begin

3288

Msg('Unable to load ' +lHdr.ImgFileName);

3289

freemem(lImg);

3290

exit;

3291

end;

3292

result := 0;

3293

for lVox := 1 to lVolVox do

3294

if (lImg^[lVox] <> 0) then

3295

inc(result);

3296

freemem(lImg);

3297

end;

3298

3299

3300

procedure TMainForm.PenalizedLogisticRegerssion1Click(Sender: TObject);

3301

label

3302

666;

3303

var

3304

lVol,lMin,lMax,lI,lFact,lnFactors,lSubj,lnSubj,lMaskVoxels,lnCrit: integer;

3305

lImageNames: TStrings;

3306

lPredictorList: TStringList;

3307

lTemp4D,lMaskname,lOutName,lStr: string;

3308

lMaskHdr: TMRIcroHdr;

3309

lMultiSymptomRA,lTempRA: singleP;

3310

//lBinomial: boolean;

3311

begin

3312

// lBinomial := false;

3313

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

3314

//next, get 1st group

3315

if not GetValX(lnSubj,lnFactors,lMultiSymptomRA,lImageNames,lnCrit{,binom},lPredictorList) then

3316

goto 666;

3317

lTemp4D := CreateDecompressed4D(lImageNames);

3318

if (lnSubj < 2) or (lnFactors < 1) then begin

3319

showmessage('This analysis requires at least 2 participants and one factor');

3320

goto 666;

3321

end;

3322

lMaskname := lImageNames[0];

3323

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

3324

showmessage('Error reading 1st image: '+lMaskname);

3325

goto 666;

3326

end;

3327

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

3328

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

3329

showmessage('File dimensions differ from mask.');

3330

goto 666;

3331

end;

3332

case MessageDlg('Do you want to add lesion volume as a regressor?', mtConfirmation,

3333

[mbYes, mbNo], 0) of { produce the message dialog box }

3334

mrYes: begin

3335

//add a new condition called lesionvolume - create a new larger array for data

3336

Msg('Computing lesion volumes...');

3337

lPredictorList.Add('LesionVolume');

3338

GetMem(lTempRA,lnSubj*lnFactors*sizeof(single));

3339

for lI := 1 to (lnSubj*lnFactors) do

3340

lTempRA^[lI] := lMultiSymptomRA^[lI];

3341

Freemem(lMultiSymptomRA);

3342

GetMem(lMultiSymptomRA,lnSubj*(lnFactors+1)*sizeof(single));

3343

for lI := 1 to (lnSubj*lnFactors) do

3344

lMultiSymptomRA^[lI] := lTempRA^[lI];

3345

Freemem(lTempRA);

3346

//now create the new factor

3347

lI := lnSubj*lnFactors;

3348

for lSubj := 1 to lnSubj do

3349

lMultiSymptomRA^[lI+lSubj] := ComputeLesionVolume(lImageNames[lSubj-1]);

3350

//ensure there is variability in this regressor

3351

lMin := round(lMultiSymptomRA^[lI+1]);

3352

lMax := round(lMultiSymptomRA^[lI+1]);

3353

for lSubj := 1 to lnSubj do begin

3354

lVol := round(lMultiSymptomRA^[lI+lSubj]);

3355

if lVol < lMin then lMin := lVol;

3356

if lVol > lMax then lMax := lVol;

3357

end;

3358

if (lMin < 0) then begin

3359

showmessage('Regression aborted: Error computing lesion volumes.');

3360

goto 666;

3361

end;

3362

if (lMin = lMax) then begin

3363

showmessage('Regression aborted: no variability in lesion volume.');

3364

goto 666;

3365

end;

3366

inc(lnFactors);

3367

end; //if user decides to include lesion volume

3368

end; //case

3369

3370

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

3371

if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file} then begin

3372

showmessage('Mask file size too small.');

3373

goto 666;

3374

end;

3375

lOutName := ExtractFileDirWithPathDelim(lMaskName)+'results';

3376

SaveHdrDlg.Filename := loutname;

3377

MsgClear;

3378

Msg(kVers);

3379

Msg('Firth Penalized regression is still beta software...');

3380

Msg('Number of participants: '+inttostr(lnSubj));

3381

Msg('Number of factors: '+inttostr(lnFactors));

3382

Msg('Threads: '+inttostr(gnCPUThreads));

3383

//next - header shows factor names

3384

lStr :='imagename';

3385

for lFact := 1 to lnFactors do

3386

lStr := lStr+','+lPredictorList[lFact-1];

3387

Msg(lStr);

3388

For lSubj := 1 to lnSubj do begin

3389

lStr :='';

3390

for lFact := 1 to lnFactors do begin

3391

lStr := lStr+','+realtostr(lMultiSymptomRA^[lSubj+ ((lFact-1)*lnSubj)],2);

3392

end;

3393

Msg (lImageNames.Strings[lSubj-1] + ' = '+lStr );

3394

end;

3395

Msg('Total voxels = '+inttostr(lMaskVoxels));

3396

Msg('Only testing voxels damaged in at least '+inttostr(lnCrit)+' individual[s]');

3397

Msg('Number of Lesion maps = '+inttostr(lnSubj));

3398

lOutName := lOutName+'.nii.gz';

3399

if not SaveHdrName ('Base Statistical Map', lOutName) then goto 666;

3400

3401

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

3402

showmessage('File dimensions differ from mask.');

3403

goto 666;

3404

end;

3405

FirthNPMAnalyze (lImageNames,lPredictorList,lMaskHdr,lnFactors,lnCrit, lMultiSymptomRA, lOutName);

3406

666:

3407

lImageNames.Free;

3408

lPredictorList.Free;

3409

DeleteDecompressed4D(lTemp4D);

3410

end;

3411

3412

function ComputeIntersection ( lAname,lBname: string; var lUnion,lIntersection,lAnotB,lBnotA: integer): boolean;

3413

label 667;

3414

var

3415

lOutName,lOutNameMod: string;

3416

lVolVox,lVolVoxA,lVox: integer;

3417

lImgA,lImgB: SingleP;

3418

3419

lMaskHdr: TMRIcroHdr;

3420

lA,lB: boolean;

3421

begin

3422

lUnion:= 0;

3423

lIntersection := 0;

3424

lAnotB := 0;

3425

lBnotA := 0;

3426

result := false;

3427

//read A

3428

if not NIFTIhdr_LoadHdr(lAname,lMaskHdr) then begin

3429

showmessage('Error reading image A - '+lAname);

3430

exit;

3431

end;

3432

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

3433

if (lVolVox < 1) then goto 667;

3434

getmem(lImgA,lVolVox*sizeof(single));

3435

if not LoadImg(lAname, lImgA, 1, lVolVox,round(lMaskHdr.NIFTIhdr.vox_offset),1,lMaskHdr.NIFTIhdr.datatype,lVolVox) then begin

3436

msg('Unable to load mask ' +lMaskHdr.ImgFileName);

3437

goto 667;

3438

end;

3439

lVolVoxA := lVolVox;

3440

//read B

3441

if not NIFTIhdr_LoadHdr(lBname,lMaskHdr) then begin

3442

showmessage('Error reading image B - '+lBname);

3443

exit;

3444

end;

3445

//fx(666,round(gOffsetRA[0]));

3446

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

3447

if (lVolVoxA <> lVolVox) or (lVolVox < 1) then goto 667;

3448

getmem(lImgB,lVolVox*sizeof(single));

3449

if not LoadImg(lBname, lImgB, 1, lVolVox,round(lMaskHdr.NIFTIhdr.vox_offset),1,lMaskHdr.NIFTIhdr.datatype,lVolVox) then begin

3450

msg('Unable to load mask ' +lMaskHdr.ImgFileName);

3451

goto 667;

3452

end;

3453

for lVox := 1 to lVolVox do begin

3454

lA := (lImgA^[lVox] <> 0);

3455

lB := (lImgB^[lVox] <> 0);

3456

if lA and lB then begin

3457

//fx(lVox,lImgA^[lVox],lImgB^[lVox]);

3458

inc(lIntersection);

3459

end;

3460

if lA or lB then

3461

inc(lUnion);

3462

if lA and not lB then

3463

inc(lAnotB);

3464

if lB and not lA then

3465

inc(lBnotA);

3466

3467

end;

3468

freemem(lImgA);

3469

freemem(lImgB);

3470

result := true;

3471

667:

3472

end;

3473

3474

procedure TMainForm.ZtoP1Click(Sender: TObject);

3475

var

3476

lAname,lBname: string; var lUnion,lIntersection,lAnotB,lBnotA: integer;

3477

begin

3478

//removed

3479

lAName := 'C:\mri\roc\p2.nii.gz';

3480

lBName := 'C:\mri\roc\RBD35.voi';

3481

if not ComputeIntersection ( lAName,lBName,lUnion,lIntersection,lAnotB,lBnotA) then

3482

Msg('Error');

3483

Msg( lAName+' '+lBName+' I'+inttostr(lIntersection)+' U'+inttostr(lUnion)+' AnotB'+inttostr(lAnotB)+' BnotA'+inttostr(lBnotA));

3484

3485

end;

3486

3487

3488

procedure TMainForm.ComputeIntersectionandUnion1Click(Sender: TObject);

3489

label

3490

666;

3491

var

3492

lUnion,lIntersection,lAnotB,lBnotA,

3493

lnSubj,lSubj,lMaskVoxels,lnAdditionalFactors: integer;

3494

lImageNames: TStrings;

3495

lMaskname,

3496

lStr,lOutName: string;

3497

lMaskHdr: TMRIcroHdr;

3498

X: PMatrix;

3499

begin

3500

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

3501

MsgClear;

3502

Msg(kVers);

3503

Msg('Compute intersection [A and B] and union [A or B] for a series of images');

3504

3505

3506

if not ReadPairedFilenamesReg(lImageNames,X,lnAdditionalFactors) then exit;

3507

lnSubj :=lImageNames.Count div 2;

3508

if lnAdditionalFactors > 1 then

3509

DelMatrix(X, lnAdditionalFactors, lnSubj);

3510

3511

lMaskname :=lImageNames[0];

3512

if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin

3513

showmessage('Error reading first image.');

3514

goto 666;

3515

end;

3516

lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);

3517

if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file} then begin

3518

showmessage('Image file size too small.');

3519

goto 666;

3520

end;

3521

3522

if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin

3523

showmessage('File dimensions differ from first image.');

3524

goto 666;

3525

end;

3526

3527

3528

Msg ('n Subjects = '+inttostr(lnSubj));

3529

for lSubj := 0 to (lnSubj-1) do begin

3530

lStr := 'A=,'+lImageNames[lSubj]+',B=,'+lImageNames[lSubj+lnSubj];

3531

ComputeIntersection ( lImageNames[lSubj],lImageNames[lSubj+lnSubj],lUnion,lIntersection,lAnotB,lBnotA);

3532

lStr := lStr + ',A and B=,'+inttostr(lIntersection);

3533

lStr := lStr + ',A or B=,'+inttostr(lUnion);

3534

lStr := lStr + ',A not B=,'+inttostr(lAnotB);

3535

lStr := lStr + ',B not A=,'+inttostr(lBnotA);

3536

3537

3538

3539

3540

Msg(lStr);

3541

end;

3542

3543

//Msg('Mask = '+lMaskname);

3544

//Msg('Total voxels = '+inttostr(lMaskVoxels));

3545

Msg('Number of observations = '+inttostr(lnSubj));

3546

666:

3547

lImageNames.Free;

3548

end; //compute intersection and union

3549

3550

3551

procedure TMainForm.ROCbinomialdeficit1Click(Sender: TObject);

3552

begin

3553

testROC;

3554

end;

3555

3556

procedure TMainForm.ROCcontinuousdeficit1Click(Sender: TObject);

3557

begin

3558

testROC2;

3559

end;

3560

3561

function isBinom ( lRA: singleP; lnObs: integer): boolean;

3562

var

3563

lI: integer;

3564

begin

3565

result := false;

3566

if lnObs < 1 then exit;

3567

for lI := 1 to lnObs do

3568

if (lRA^[lI] <> 0) and (lRA^[lI] <> 1) then

3569

exit;

3570

result := true;

3571

end;

3572

3573

procedure Means ( lBinomRA,lContRA: singleP; lnObs: integer);

3574

var

3575

lI,ln0: integer;

3576

lMeans0, lMeans1: double;

3577

begin

3578

lMeans0 := 0;

3579

lMeans1 := 0;

3580

ln0 := 0;

3581

if lnObs < 1 then exit;

3582

for lI := 1 to lnObs do begin

3583

if (lBinomRA^[lI] = 0) then begin

3584

inc(ln0);

3585

lMeans0 := lMeans0 + lContRA^[lI];

3586

end else

3587

lMeans1 := lMeans1 + lContRA^[lI];

3588

end;

3589

if ln0 > 0 then

3590

lMeans0 := lMeans0 / ln0;

3591

if ln0 < lnObs then

3592

lMeans1 := lMeans1 / (lnObs-ln0);

3593

npmform.MainForm.memo1.lines.add('mean volume for '+inttostr(ln0)+' people who scored 0 is = '+floattostr(lmeans0));

3594

npmform.MainForm.memo1.lines.add('mean volume for '+inttostr(lnObs-ln0)+' people who scored 1 is = '+floattostr(lmeans1));

3595

3596

end;

3597

3598

function AUCbinomcontT (lBinomdataRA,lContdataRA: singlep; lnSubj :integer; var lT: double): double;

3599

var

3600

lIn : DoubleP0;

3601

lnGroup0,lnGroup1,lI: integer;

3602

begin

3603

result := 0.5;

3604

if lnSubj < 1 then

3605

exit;

3606

Getmem(lIn,lnSubj*sizeof(double));

3607

lnGroup0 := 0;

3608

lnGroup1 := 0;

3609

for lI := 1 to lnSubj do begin

3610

if lBinomdataRA^[lI] = 0 then begin

3611

lIn^[lnGroup0] := lContdataRA^[lI];

3612

inc (lnGroup0);

3613

end else begin

3614

inc (lnGroup1);

3615

lIn^[lnSubj-lnGroup1] := lContdataRA^[lI];

3616

3617

end;

3618

end;

3619

result := continROC (lnSubj, lnGroup0, lIn);

3620

TStat2 (lnSubj, lnGroup0, lIn,lT);

3621

freemem(lIn);

3622

end;

3623

3624

3625

procedure Contrast(lBehavName,lROIname: string; lBehavRA,lLesionVolRA: singleP; lnSubj: integer);

3626

var

3627

lDF: integer;

3628

lROC,lT,lP: double;

3629

begin

3630

if isBinom (lBehavRA,lnSubj) then begin

3631

lROC := AUCbinomcontT (lBehavRA,lLesionVolRA, lnSubj,lT);

3632

lDF := lnSubj-2;

3633

lP := pTdistr(lDF,lT);

3634

Means ( lBehavRA,lLesionVolRA, lnSubj);

3635

3636

npmform.MainForm.memo1.lines.add('ROI=,'+lROIname+',Behav=,'+lBehavName+', Area Under Curve=,'+floattostr(lROC)+', T('+inttostr(lDF)+')=,'+floattostr(lT)+',p<,'+floattostr(lp));

3637

end else begin

3638

lROC := AUCcontcont (lBehavRA,lLesionVolRA, lnSubj);

3639

npmform.MainForm.memo1.lines.add('ROI=,'+lROIname+',Behav=,'+lBehavName+', Area Under Curve = '+floattostr(lROC));

3640

end;

3641

//xxx

3642

end;

3643

3644

function ComputeOverlap ( lROIname: string; var lLesionNames: TStrings; var lROIvol: double; lFracROIinjured: singlep): boolean;

3645

label 667;

3646

var

3647

lName: string;

3648

lSum: double;

3649

lLesion,lnLesions,lVolVox,lVolVoxA,lVox: integer;

3650

lROIImg,lImgB: SingleP;

3651

lMaskHdr: TMRIcroHdr;

3652

begin

3653

lROIvol := 0;

3654

result := false;

3655

lnLesions := lLesionNames.count;

3656

if lnLesions < 1 then begin

3657

showmessage('Error: no lesion names');

3658

exit;

3659

end;

3660

for lLesion := 1 to lnLesions do

3661

lFracROIinjured^[lLesion] := 0;

3662

//read A

3663

if not NIFTIhdr_LoadHdr(lROIname,lMaskHdr) then begin

3664

showmessage('Error reading ROI - '+lROIname);

3665

exit;

3666

end;

3667

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

3668

if (lVolVox < 1) then begin

3669

showmessage('Error with Mask voxels ' + inttostr(lVolVox));

3670

exit;

3671

end;

3672

getmem(lROIImg,lVolVox*sizeof(single));

3673

getmem(lImgB,lVolVox*sizeof(single));

3674

if not LoadImg(lROIname, lROIImg, 1, lVolVox,round(lMaskHdr.NIFTIhdr.vox_offset),1,lMaskHdr.NIFTIhdr.datatype,lVolVox) then begin

3675

MainForm.NPMmsg('Unable to load mask ' +lMaskHdr.ImgFileName);

3676

goto 667;

3677

end;

3678

lVolVoxA := lVolVox;

3679

for lVox := 1 to lVolVox do

3680

if (lROIImg^[lVox] > 0) then

3681

lROIvol := lROIvol +lROIImg^[lVox];

3682

//read Lesion

3683

if lROIvol < 1 then begin

3684

MainForm.NPMmsg('ROI volume < 1');

3685

goto 667;

3686

end;

3687

//for each lesion

3688

//MainForm.NPMmsg('Compute overlap started '+inttostr(lnLesions)+' '+floattostr(lROIvol)+' '+inttostr(lVolVoxA));

3689

MainForm.ProgressBar1.Position := 0;

3690

for lLesion := 1 to lnLesions do begin

3691

MainForm.ProgressBar1.Position := round((lLesion/lnLesions)*100) ;

3692

lSum := 0;

3693

lName := lLesionNames.Strings[lLesion-1];

3694

if not NIFTIhdr_LoadHdr(lName,lMaskHdr) then begin

3695

showmessage('Error reading lesion - '+lName);

3696

exit;

3697

end;

3698

lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];

3699

if (lVolVoxA <> lVolVox) or (lVolVox < 1) then begin

3700

MainForm.NPMmsg('Volume does not have expected number of voxels ' +lMaskHdr.ImgFileName +'<>'+lROIname+ ' found ' +inttostr(lVolVox)+' expected '+inttostr(lVolVoxA));

3701

goto 667;

3702

end;

3703

if not LoadImg(lName, lImgB, 1, lVolVox,round(lMaskHdr.NIFTIhdr.vox_offset),1,lMaskHdr.NIFTIhdr.datatype,lVolVox) then begin

3704

MainForm.NPMmsg('Unable to load mask ' +lMaskHdr.ImgFileName);

3705

goto 667;

3706

end;

3707

for lVox := 1 to lVolVox do begin

3708

//if {(lImgB^[lVox] <> 0) and} (lROIImg^[lVox] <> 0) then fx(lROIImg^[lVox]);

3709

3710

if (lROIImg^[lVox] > 0) and (lImgB^[lVox] <> 0) then

3711

lSum := lSum + lROIImg^[lVox];

3712

end;

3713

lFracROIinjured^[lLesion] := lSum/lROIvol;

3714

end;//for each lesion

3715

result := true;

3716

MainForm.ProgressBar1.Position := 0;

3717

3718

(*for lLesion := 1 to lnLesions do begin

3719

if lFracROIinjured^[lLesion] > 0 then

3720

fx( lFracROIinjured^[lLesion], lLesion);

3721

end; *)

3722

3723

667:

3724

3725

freemem(lImgB);

3726

freemem(lROIImg);

3727

end;

3728

3729

3730

procedure TMainForm.ROIanalysis1Click(Sender: TObject);

3731

label

3732

666;

3733

var

3734

lROI,lnROI,lVol,lMin,lMax,lI,lFact,lnFactors,lSubj,lnSubj,lMaskVoxels,lnCrit: integer;

3735

lROInames,lImageNames: TStrings;

3736

lPredictorList: TStringList;

3737

lTemp4D,lOutName,lStr: string;

3738

lBehav: single;

3739

lROIvolRA: doubleP;

3740

lMultiSymptomRA,lLesionVolRA,lBehavRA: singleP;

3741

lError: boolean;

3742

begin

3743

//OpenHdrDlg.Files.Add( 'C:\lesions\vois\frontal2.voi');

3744

//OpenHdrDlg.Files.Add( 'C:\lesions\vois\posterior2.voi');

3745

//MainForm.OpenHdrDlg.InitialDir := 'C:\lesions';

3746

if not OpenDialogExecute('Select regions of interest',true,false,kImgPlusVOIFilter) then begin

3747

showmessage('NPM aborted: file selection failed.');

3748

exit;

3749

end; //if not selected

3750

lROInames:= TStringList.Create;

3751

lROInames.addstrings(OpenHdrDlg.Files);

3752

lnROI := lROINames.Count;

3753

if lnROI < 1 then begin

3754

showmessage('You need to select at least one ROI.');

3755

exit;

3756

end;

3757

lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???

3758

if not GetValX(lnSubj,lnFactors,lMultiSymptomRA,lImageNames,lnCrit,lPredictorList) then

3759

goto 666;

3760

lTemp4D := CreateDecompressed4D(lImageNames);

3761

if (lnSubj < 1) or (lnFactors < 1) then begin

3762

showmessage('This analysis requires at least 1 participant and one factor');

3763

goto 666;

3764

end;

3765

MsgClear;

3766

Msg(kVers);

3767

MainForm.NPMmsg('Analysis began = ' +TimeToStr(Now));

3768

Msg('Number of participants: '+inttostr(lnSubj));

3769

Msg('Number of factors: '+inttostr(lnFactors));

3770

Msg('Number of Lesion maps = '+inttostr(lnSubj));

3771

//next - header shows factor names

3772

lStr :='imagename';

3773

for lFact := 1 to lnFactors do

3774

lStr := lStr+','+lPredictorList[lFact-1];

3775

for lROI := 1 to lnROI do

3776

lStr := lStr+','+lROInames[lROI-1];

3777

Msg(lStr);

3778

lError := false;

3779

Getmem(lROIVolRA, lnSubj*lnROI*sizeof(double));

3780

Getmem(lLesionVolRA, lnSubj*lnROI*sizeof(single));

3781

Getmem(lBehavRA, lnSubj*lnFactors*sizeof(single));

3782

for lROI := 1 to lnROI do begin

3783

//if not ComputeIntersection ( lImageNames.Strings[lSubj-1],lROInames[lROI-1],lUnion,lIntersection,lAnotB,lBnotA) then

3784

if not ComputeOverlap (lROInames[lROI-1],lImageNames, lROIvolRA^[lROI], singlep(@lLesionVolRA^[((lROI-1)*lnSubj)+1])) then begin

3785

MainForm.NPMmsg('Error computing overlap');

3786

goto 666;

3787

end;

3788

end;

3789

For lSubj := 1 to lnSubj do begin

3790

lStr :='';

3791

for lFact := 1 to lnFactors do begin

3792

lBehav := lMultiSymptomRA^[lSubj+ ((lFact-1)*lnSubj)];

3793

lStr := lStr+','+realtostr(lBehav,2);

3794

lBehavRA^[((lFact-1)*lnSubj) +lSubj] := lBehav;

3795

end;

3796

for lROI := 1 to lnROI do

3797

lStr := lStr+','+floattostr(lLesionVolRA^[((lROI-1)*lnSubj) +lSubj]);

3798

Msg (lImageNames.Strings[lSubj-1] + ' = '+lStr );

3799

end;

3800

(* For lSubj := 1 to lnSubj do begin

3801

lStr :='';

3802

for lFact := 1 to lnFactors do begin

3803

lBehav := lMultiSymptomRA^[lSubj+ ((lFact-1)*lnSubj)];

3804

lStr := lStr+','+realtostr(lBehav,2);

3805

lBehavRA^[((lFact-1)*lnSubj) +lSubj] := lBehav;

3806

end;

3807

for lROI := 1 to lnROI do begin

3808

if ComputeIntersection ( lImageNames.Strings[lSubj-1],lROInames[lROI-1],lUnion,lIntersection,lAnotB,lBnotA) then begin

3809

lStr := lStr+','+inttostr(lIntersection);

3810

lLesionVolRA^[((lROI-1)*lnSubj) +lSubj] := lIntersection;

3811

end else begin

3812

lError:= true;

3813

lStr := lStr+',error';

3814

end;

3815

//Msg( lImageNames.Strings[lSubj-1]+' '+lROInames[lROI-1]+' I'+inttostr(lIntersection)+' U'+inttostr(lUnion)+' AnotB'+inttostr(lAnotB)+' BnotA'+inttostr(lBnotA));

3816

3817

end;

3818

Msg (lImageNames.Strings[lSubj-1] + ' = '+lStr );

3819

end;*)

3820

for lROI := 1 to lnROI do begin

3821

for lFact := 1 to lnFactors do begin

3822

Contrast(lPredictorList[lFact-1],lROInames[lROI-1],singlep(@lBehavRA^[((lFact-1)*lnSubj)+1]),singlep(@lLesionVolRA^[((lROI-1)*lnSubj)+1]),lnSubj);//,((lFact-1)*lnSubj),((lROI-1)*lnSubj));

3823

end; //for each factor

3824

end; //for each ROI

3825

for lROI := 1 to lnROI do begin

3826

Msg( lROInames[lROI-1] +' volume = '+floattostr(lROIvolRA^[lROI]) )

3827

end; //for each ROI

3828

Freemem(lLesionVolRA);

3829

Freemem(lBehavRA);

3830

Freemem(lROIvolRA);

3831

3832

666:

3833

lROInames.free;

3834

lImageNames.Free;

3835

lPredictorList.Free;

3836

DeleteDecompressed4D(lTemp4D);

3837

MainForm.NPMmsg('Analysis finished = ' +TimeToStr(Now));

3838

end;

3839

3840

3841

procedure TMainForm.Masked1Click(Sender: TObject);

3842

var

3843

lFilename,lMaskname: string;

3844

lPos: Integer;

3845

begin

3846

MsgClear;

3847

Msg(GetKVers);

3848

if not OpenDialogExecute('Select brain mask ',false,false,kImgFilter) then begin

3849

showmessage('NPM aborted: mask selection failed.');

3850

exit;

3851

end; //if not selected

3852

lMaskname := OpenHdrDlg.Filename;

3853

if not OpenDialogExecute('Select images for intensity normalization',true,false,kImgFilter) then begin

3854

showmessage('NPM aborted: file selection failed.');

3855

exit;

3856

end; //if not selected

3857

if OpenHdrDlg.Files.Count < 1 then

3858

exit;

3859

3860

for lPos := 1 to OpenHdrDlg.Files.Count do begin

3861

lFilename := OpenHdrDlg.Files[lPos-1];

3862

balance(lFilename,lMaskname,(Sender as TMenuItem).tag);

3863

end;

3864

end;

3865

3866

function Binarize (var lImageName:String; lNonZeroVal: integer; lZeroThresh: boolean): boolean;

3867

var

3868

lImg8: ByteP;

3869

lImg: SingleP;

3870

lHdr: TMRIcroHdr;

3871

lVolVox,lVox: integer;

3872

lMin,lMax: single;

3873

lModeLo,lModeHi,lIntercept,lSlope: single;

3874

lOutNameMod: string;

3875

begin

3876

//lOutName := lMaskHdr.ImgFileName;

3877

result := false;

3878

//if not SaveHdrName ('Statistical Map', lOutName) then exit;

3879

if not NIFTIhdr_LoadHdr(lImageName,lHdr) then begin

3880

showmessage('Error reading '+lImageName);

3881

exit;

3882

end;

3883

lVolVox := lHdr.NIFTIhdr.dim[1]*lHdr.NIFTIhdr.dim[2]* lHdr.NIFTIhdr.dim[3];

3884

if (lVolVox < 1) then exit;

3885

getmem(lImg,lVolVox*sizeof(single));

3886

getmem(lImg8,lVolVox*sizeof(byte));

3887

if not LoadImg(lHdr.ImgFileName, lImg, 1, lVolVox,round(lHdr.NIFTIhdr.vox_offset),1,lHdr.NIFTIhdr.datatype,lVolVox) then begin

3888

Msg('Unable to load ' +lHdr.ImgFileName);

3889

exit;

3890

end;

3891

lHdr.NIFTIhdr.scl_slope := 1;

3892

lHdr.NIFTIhdr.scl_inter := 0;

3893

if lZeroThresh then begin

3894

lOutNameMod := ChangeFilePrefixExt(lImageName,'i','.nii');

3895

3896

lMin := 0;

3897

lMax := 0

3898

end else begin

3899

lOutNameMod := ChangeFilePrefixExt(lImageName,'i','.voi');

3900

3901

lMin := lIMg^[1];

3902

for lVox := 1 to lVolVox do

3903

if lImg^[lVox] < lMin then lMin := lIMg^[lVox];

3904

3905

lMax := lIMg^[1];

3906

for lVox := 1 to lVolVox do

3907

if lImg^[lVox] > lMax then lMax := lIMg^[lVox];

3908

for lVox := 1 to lVolVox do

3909

lImg8^[lVox] := 0;

3910

lMax := ((lMax-lMin) / 2)+lMin;

3911

end;

3912

for lVox := 1 to lVolVox do

3913

if lImg^[lVox] > lMax then

3914

lImg8^[lVox] := lNonZeroVal;

3915

Msg('Creating ' +lOutNameMod+' Threshold = '+floattostr(lMax));

3916

NIFTIhdr_SaveHdrImg8(lOutNameMod,lHdr.NIFTIhdr,true,not IsNifTiMagic(lHdr.NIFTIhdr),true,lImg8,1);

3917

freemem(lIMg8);

3918

freemem(lImg);

3919

3920

end;

3921

3922

3923

procedure TMainForm.Binarizeimages1Click(Sender: TObject);

3924

var

3925

lFilename: string;

3926

lPos: Integer;

3927

begin

3928

MsgClear;

3929

Msg(GetKVers);

3930

if not OpenDialogExecute('Select images for intensity normalization',true,false,kImgFilter) then begin

3931

showmessage('NPM aborted: file selection failed.');

3932

exit;

3933

end; //if not selected

3934

if OpenHdrDlg.Files.Count < 1 then

3935

exit;

3936

for lPos := 1 to OpenHdrDlg.Files.Count do begin

3937

lFilename := OpenHdrDlg.Files[lPos-1];

3938

Binarize(lFilename,1,false);

3939

//Binarize (var lImageName:String; lNonZeroVal: integer; lZeroThresh: boolean): boolean;

3940

end;

3941

Msg('Done');

3942

end;

3943

3944

procedure TMainForm.Resliceimagetoneworientationandboundingbox1Click(

3945

Sender: TObject);

3946

begin

3947

(* var

3948

lSourcename,lTargetName: string;

3949

lPos: integer;

3950

begin

3951

MsgClear;

3952

Msg(GetKVers);

3953

Msg('This function will transform a source image to match a target image.');

3954

Msg(' The resliced image will have the voxel size, orientation and bounding box of the target.');

3955

Msg('You will be prompted to select a target image as well as source images.');

3956

Msg(' The resliced image will have the prefix ''r'', e.g. c:\source.nii -> c:\rsource.nii.');

3957

if not OpenDialogExecute('Select target image (source images will be resliced to match target)',false,false,kImgFilter) then begin

3958

showmessage('reslice aborted: target selection failed.');

3959

exit;

3960

end; //if not selected

3961

lTargetName := OpenHdrDlg.Filename;

3962

if not OpenDialogExecute('Select source images for reslicing',true,false,kImgFilter) then begin

3963

showmessage('reslice aborted: source selection failed.');

3964

exit;

3965

end; //if not selected

3966

if OpenHdrDlg.Files.Count < 1 then

3967

exit;

3968

for lPos := 1 to OpenHdrDlg.Files.Count do begin

3969

lSourcename := OpenHdrDlg.Files[lPos-1];

3970

3971

xxBinarize(lFilename);

3972

end;

3973

Msg('Done');

3974

3975

end;

3976

3977

3978

procedure TMainForm.Setnonseroto1001Click(Sender: TObject);

3979

var

3980

lFilename: string;

3981

lPos: Integer;

3982

begin

3983

MsgClear;

3984

Msg(GetKVers);

3985

if not OpenDialogExecute('Select images for intensity normalization',true,false,kImgFilter) then begin

3986

showmessage('NPM aborted: file selection failed.');

3987

exit;

3988

end; //if not selected

3989

if OpenHdrDlg.Files.Count < 1 then

3990

exit;

3991

for lPos := 1 to OpenHdrDlg.Files.Count do begin

3992

lFilename := OpenHdrDlg.Files[lPos-1];

3993

Binarize(lFilename,100,true);

3994

//Binarize (var lImageName:String; lNonZeroVal: integer; lZeroThresh: boolean): boolean;

3995

end;

3996

end;

3997

3998

procedure TMainForm.Savetext1Click(Sender: TObject);

3999

begin

4000

SaveHdrDlg.Title := 'Save file as comma separated values (to open with Excel)';

4001

SaveHdrDlg.Filter := 'Comma Separated (*.csv)|*.csv|Text (*.txt)|*.txt';

4002

SaveHdrDlg.DefaultExt := '*.csv';

4003

if not SaveHdrDlg.Execute then exit;

4004

Memo1.Lines.SaveToFile(SaveHdrDlg.Filename);

4005

end;

4006

4007

procedure TMainForm.AnaCOMmenuClick(Sender: TObject);

4008

begin

4009

{$IFDEF compileANACOM}

4010

DoAnaCOM;

4011

{$ENDIF}

4012

end;

4013

4014

procedure TMainForm.MonteCarloSimulation1Click(Sender: TObject);

4015

begin

4016

{$IFDEF benchmark}

4017

LesionMonteCarlo (false,true,true);

4018

{$ENDIF}

4019

end;

4020

4021

function TMainForm.MakeSubtract (lPosName,lNegName: string): boolean;

4022

var

4023

lNegImg,lImg,lOutImg: SingleP;

4024

lHdr,lNegHdr: TMRIcroHdr;

4025

lVolVox,lVox: integer;

4026

lOutNameMod: string;

4027

begin

4028

result := false;

4029

if not NIFTIhdr_LoadHdr(lPosName,lHdr) then begin

4030

showmessage('Error reading '+lPosName);

4031

exit;

4032

end;

4033

lVolVox := lHdr.NIFTIhdr.dim[1]*lHdr.NIFTIhdr.dim[2]* lHdr.NIFTIhdr.dim[3];

4034

if (lVolVox < 1) then exit;

4035

getmem(lImg,lVolVox*sizeof(single));

4036

if not LoadImg(lHdr.ImgFileName, lImg, 1, lVolVox,round(lHdr.NIFTIhdr.vox_offset),1,lHdr.NIFTIhdr.datatype,lVolVox) then begin

4037

Msg('Unable to load ' +lHdr.ImgFileName);

4038

exit;

4039

end;

4040

4041

if not NIFTIhdr_LoadHdr(lNegName,lNegHdr) then begin

4042

showmessage('Error reading '+lNegName);

4043

exit;

4044

end;

4045

if lVolVox <> (lNegHdr.NIFTIhdr.dim[1]*lNegHdr.NIFTIhdr.dim[2]* lNegHdr.NIFTIhdr.dim[3]) then begin

4046

showmessage('Volumes differ');

4047

exit;

4048

4049

end;

4050

getmem(lImg,lVolVox*sizeof(single));

4051

if not LoadImg(lHdr.ImgFileName, lImg, 1, lVolVox,round(lHdr.NIFTIhdr.vox_offset),1,lHdr.NIFTIhdr.datatype,lVolVox) then begin

4052

Msg('Unable to load ' +lHdr.ImgFileName);

4053

exit;

4054

end;

4055

getmem(lNegImg,lVolVox*sizeof(single));

4056

if not LoadImg(lNegHdr.ImgFileName, lNegImg, 1, lVolVox,round(lNegHdr.NIFTIhdr.vox_offset),1,lNegHdr.NIFTIhdr.datatype,lVolVox) then begin

4057

Msg('Unable to load ' +lNegHdr.ImgFileName);

4058

exit;

4059

end;

4060

getmem(lOutImg,lVolVox*sizeof(single));

4061

for lVox := 1 to lVolVox do

4062

lOutImg^[lVox] := lImg^[lVox] - lNegImg^[lVox];

4063

4064

4065

lHdr.NIFTIhdr.scl_slope := 1;

4066

lHdr.NIFTIhdr.scl_inter := 0;

4067

lOutNameMod := ChangeFilePrefixExt(lPosName,'subtract_','.hdr');

4068

Msg(lPosName+' - ' + lNegName+ ' = '+lOutNameMod);

4069

NIFTIhdr_SaveHdrImg(lOutNameMod,lHdr.NIFTIhdr,true,not IsNifTiMagic(lHdr.NIFTIhdr),true,lOutImg,1);

4070

//end optional

4071

//NIFTIhdr_SaveHdr(lHdr.HdrFilename,lHdr.NIFTIhdr,true,not IsNifTiMagic(lHdr.NIFTIhdr));

4072

4073

freemem(lImg);

4074

freemem(lOutImg);

4075

freemem(lNegImg);

4076

end;//makesubtract

4077

4078

procedure TMainForm.Subtract1Click(Sender: TObject);

4079

var

4080

lPosName,lNegName: string;

4081

begin

4082

if not OpenDialogExecute('Select positive',false,false,kImgPlusVOIFilter) then

4083

exit;

4084

lPosName := OpenHdrDlg.FileName;

4085

if not OpenDialogExecute('Select negative',false,false,kImgPlusVOIFilter) then

4086

exit;

4087

lNegName := OpenHdrDlg.FileName;

4088

MakeSubtract (lPosName,lNegName);

4089

4090

end;

4091

4092

4093

4094

4095

4096

4097

procedure TMainForm.LogPtoZ1Click(Sender: TObject);

4098

var

4099

lFilename: string;

4100

lPos: Integer;

4101

begin

4102

MsgClear;

4103

Msg(GetKVers);

4104

if not OpenDialogExecute('Select images for intensity normalization',true,false,kImgFilter) then begin

4105

showmessage('NPM aborted: file selection failed.');

4106

exit;

4107

end; //if not selected

4108

if OpenHdrDlg.Files.Count < 1 then

4109

exit;

4110

for lPos := 1 to OpenHdrDlg.Files.Count do begin

4111

lFilename := OpenHdrDlg.Files[lPos-1];

4112

//LogPToZ(lFilename,1,false);

4113

4114

end;

4115

Msg('Done');

4116

end;

4117

4118

{$IFDEF UNIX}

4119

initialization

4120

{$I npmform.lrs}

4121

{$ELSE} //not unix: windows

4122

initialization

4123

{$IFDEF FPC}

4124

{$I npmform.lrs}

4125

{$ENDIF}//FPC

4126

OleInitialize(nil);

4127

4128

finalization

4129

OleUninitialize

4130

{$ENDIF} //Windows

4131

4132

end.

4133

Older »