File: turbolesion.pas

package info (click to toggle)
mricron 0.20140804.1~dfsg.1-1
  • links: PTS, VCS
  • area: main
  • in suites: jessie, jessie-kfreebsd
  • size: 13,480 kB
  • ctags: 8,011
  • sloc: pascal: 114,853; sh: 49; makefile: 32
file content (491 lines) | stat: -rwxr-xr-x 21,735 bytes parent folder | download | duplicates (7)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
unit turbolesion;
interface
{$H+}                                                                                                                                 
uses
  define_types,SysUtils,
part,StatThds,statcr,StatThdsUtil,Brunner,DISTR,nifti_img, hdr,
   Messages,  Classes, Graphics, Controls, Forms, Dialogs,
StdCtrls,  ComCtrls,ExtCtrls,Menus, overlap,ReadInt,lesion_pattern,stats,LesionStatThds,nifti_hdr,

{$IFDEF FPC} LResources,gzio2,
{$ELSE} gziod,associate,{$ENDIF}   //must be in search path, e.g. C:\pas\mricron\npm\math
{$IFNDEF UNIX} Windows, {$ENDIF}
upower,firthThds,firth,IniFiles,cpucount,userdir,math,
regmult,utypes;
Type
  TLDMPrefs = record
         NULP,BMtest,Ttest,Ltest: boolean;
         CritPct,nCrit,nPermute,Run: integer;
         ValFilename, OutName, ExplicitMaskName: string;
  end;
function TurboLDM (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; var lPrefs: TLDMPrefs ; var lSymptomRA: SingleP;var lFactname,lOutName: string): boolean;



implementation

uses npmform;

(*procedure Debog (var lSumImg: Smallintp; lVox: integer);
var
   lInName : string;
   lFData: file;
begin
         lInName := 'c:\16.img';
	 assignfile(lFdata,lInName);
	 filemode := 2;
	 Rewrite(lFdata,lVox*sizeof(smallint));
	 BlockWrite(lFdata,lSumImg^, 1  {, NumWritten});
	 closefile(lFdata);
end;*)

function MakeSum (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; var lSumImg: Smallintp): boolean;
//if successful, you MUST freemem(lSumImg)...
label
	667;
var
	lVolVox,lVox,lImg,lPosPct: integer;
        lVolImg: byteP;

begin
        result := false;
	lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];
	if (lVolVox < 1) then exit;
        getmem(lVolImg,lVolVox* sizeof(byte));
        getmem(lSumImg,lVolVox* sizeof(smallint));
        for lVox := 1 to lVolVox do //June 2009 init array
               lSumImg^[lVox] := 0;
(*        for lVox := 1 to lVolVox do
            if lVolImg^[lVox] <> 0 then
               lSumImg^[lVox] := lSumImg^[lVox]+1;*)
        for lImg := 1 to lImages.Count do begin
                lPosPct := round(100*(lImg / lImages.Count));
                MainForm.ProgressBar1.Position := lPosPct;
                Application.Processmessages;
                if not LoadImg8(lImages[lImg-1], lVolImg, 1, lVolVox,round(gOffsetRA[lImg]),1,gDataTypeRA[lImg],lVolVox) then
                   goto 667;
                for lVox := 1 to lVolVox do
                    if lVolImg^[lVox] <> 0 then
                       lSumImg^[lVox] := lSumImg^[lVox]+1;
	end;//for each image
	MainForm.NPMmsg('Sum image finished = ' +TimeToStr(Now));
        MainForm.ProgressBar1.Position := 0;
        //Debog(lSumImg, lVolVox);
        freemem(lVolImg);
        result := true;
	exit;
667: //you only get here if you aborted ... free memory and report error
           freemem(lVolImg);
           freemem(lSumImg);
	MainForm.NPMMsg('Unable to complete analysis.');
        MainForm.ProgressBar1.Position := 0;
end;


function ThreshSumImg (var lSumImg: Smallintp; lVolVox,lThresh: integer): integer;
//sets all voxels with values < lThresh to zero, returns number of voxels to survive threshold.
var
   lPos: integer;
begin
     result := 0;
     if lVolVox < 1 then
        exit;
     for lPos := 1 to lVolVox do
         if lSumImg^[lPos] < lThresh then
            lSumImg^[lPos] := 0
         else
             inc(result);
end;

function ExplicitMaskSumImg (lMaskName: string; var lSumImg: Smallintp; lVolVox: integer): integer;
//Any voxels in MaskImg that are 0 are zeroed in the SumImg
var
   lOK: boolean;
   lPos: integer;
   lMaskHdr: TMRIcroHdr;
   lMaskData: bytep;
label
  666;
begin
     result := 0;
     if (lVolVox < 1) or (not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr)) then begin
	      MainForm.NPMmsg('Error: unable to load explicit mask named '+lMaskName);
	      exit;
     end;
     if lVolVox <> (lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3]) then begin
	      MainForm.NPMmsg('Error: data and explicit mask have different sizes '+lMaskName);
	      exit;
     end;

     getmem(lMaskData,lVolVox* sizeof(byte));
     lOK := LoadImg8(lMaskName, lMaskData, 1, lVolVox,round(lMaskHdr.NIFTIhdr.vox_offset),1,lMaskHdr.NIFTIhdr.DataType,lVolVox);
     if not lOK then goto 666;
     if lVolVox < 1 then
        exit;
     for lPos := 1 to lVolVox do
         if lMaskData^[lPos] < 1 then
            lSumImg^[lPos] := 0
         else
             inc(result);

     666:
     freemem(lMaskData);
end;

function LoadImg8Masked(lInName: string; lImgData: bytep; lMaskData: SmallIntP; lStartMaskPos, lEndMaskPos,linvox_offset,lRApos,lDataType,lVolVox: integer): boolean;
label
     111;
var
   lFullImgData: bytep;
   lMaskPos,lPos: integer;
begin
     result := false;
     if (lVolVox < 1) or (lEndMaskPos < lStartMaskPos) then
        exit;
     getmem(lFullImgData,lVolVox* sizeof(byte));
     result := LoadImg8(lInName, lFullImgData, 1, lVolVox,linvox_offset,1,lDataType,lVolVox);
     if result then begin
        lMaskPos := 0;
        for lPos := 1 to lVolVox do begin
            if lMaskData^[lPos] <> 0 then begin
                inc(lMaskPos);
                if (lMaskPos >=lStartMaskPos) then
                   lImgData^[lRApos+lMaskPos-1] := lFullImgData^[lPos];
                if lMaskPos = lEndMaskPos then goto 111;

            end;//voxel in mask
        end; //for each voxel in image

     end;//if LoadImg8 success
111:
     freemem(lFullImgData);
end;

function reformat(var lStatImg: singlep; lMaskImg: smallintp; lVolVox: integer): boolean;
var
   lPos,lStatPos,lMaskItems: integer;
begin
     result := false;
     if lVolVox < 1 then
        exit;
     lMaskItems := 0;
     for lPos := 1 to lVolVox do
         if lMaskImg^[lPos] <> 0 then
            inc(lMaskItems);
     result := true;
     if (lMaskItems < 1) or (lMaskItems >= lVolVox) then
        exit;//no need to reformat
     //note that we do this in descending order, so we do not overwrite...
     lStatPos := lMaskItems;
     for lPos := lVolVox downto 1 do
         if lMaskImg^[lPos] <> 0 then begin
            lStatImg^[lPos] := lStatImg^[lStatPos];
            dec(lStatPos);
         end else
             lStatImg^[lPos] := 0;
end;//reformat


function NULPcount (lPlankImg: bytep; lVoxPerPlank,lImagesCount: integer; var lUniqueOrders: integer; var lOverlapRA: Overlapp): boolean;
procedure CheckOrder(var lObservedOrder: TLesionPattern);
var
   lInc: integer;
begin
     if lUniqueOrders > 0 then begin //see if this is unique
          for lInc := 1 to lUniqueOrders do
              if SameOrder(lObservedOrder,lOverlapRA^[lInc],lImagesCount) then
                 exit; //not unique
     end; //UniqueOrders > 0
     //if we have not exited yet, we have found a new ordering!
     lUniqueOrders := lUniqueOrders + 1;
     lOverlapRA^[lUniqueOrders] := lObservedOrder;
end;

var
   lVox,lPlankImgPos,lPos: integer;
   lOrder,lPrevOrder: TLesionPattern;
begin
     result := false;
     lPrevOrder := EmptyOrder;//impossible: forces first voxel of each order to be checked
     for lVox := 1 to lVoxPerPlank do begin
         (*if (lVox mod lVoxPerPlankDiv10) = 0 then begin
                       MainForm.ProgressBar1.Position := (lVox div lVoxPerPlankDiv10)*10;
                       MainForm.Refresh;
                       Application.processmessages;
                    end;*)
         lOrder := EmptyOrder;
         lPlankImgPos := 0;
         //lnDeficits := 0;
         for lPos := 1 to lImagesCount do begin
             if (lPlankImg^[lPlankImgPos + lVox] > 0) then begin
                //inc(lnDeficits);
                SetBit(lPos,lOrder);
             end;
             lPlankImgPos := lPlankImgPos + lVoxPerPlank;
         end;
         //if  (lnDeficits >= lminDeficits) then begin //this is different from the last voxel: perhaps this is a new ordering
             if (not SameOrder(lOrder,lPrevOrder,lImagesCount)) then
                CheckOrder(lOrder);
             //inc(lnVoxels);
         //end;//nDeficies
         lPrevOrder := lOrder;
     end;//for lVox
     result := true;
end;

procedure PtoZpermute (lnPermute: integer; lPermuteMaxT, lPermuteMinT: singlep);
var
   lPos: integer;
   lVal : single;
begin
     if lPos < 1 then exit;
     for lPos := 1 to lnPermute do begin
            if (lPermuteMinT^[lPos] > 1.1) or (lPermuteMinT^[lPos] < -1.1) then
               lPermuteMinT^[lPos] := 0.5;
            if (lPermuteMaxT^[lPos] > 1.1) or (lPermuteMaxT^[lPos] < -1.1) then
               lPermuteMaxT^[lPos] := 0.5;
            lVal := lPermuteMaxT^[lPos];
            lPermuteMaxT^[lPos] := lPermuteMinT^[lPos];
            lPermuteMinT^[lPos] := lVal;
            if lPermuteMaxT^[lPos] < 0 then
			lPermuteMaxT^[lPos] := -pNormalInv(abs(lPermuteMaxT^[lPos]))
            else
			lPermuteMaxT^[lPos] := pNormalInv(lPermuteMaxT^[lPos]);
            if lPermuteMinT^[lPos] < 0 then
			lPermuteMinT^[lPos] := -pNormalInv(abs(lPermuteMinT^[lPos]))
            else
			lPermuteMinT^[lPos] := pNormalInv(lPermuteMinT^[lPos]);
     end;
end;


function TurboLDM (var lImages: TStrings; var lMaskHdr: TMRIcroHdr;var lPrefs: TLDMPrefs ; var lSymptomRA: SingleP;var lFactname,lOutName: string): boolean;
label
	123,667;
var
	lOutNameMod: string;
	lStatHdr: TNIfTIhdr;
        lThreshFDR,lThreshPermute,lThreshBonf,lThreshNULP :double;
	lObsp: pointer;
	lObs: Doublep0;
        lRanOrderp: pointer;
        lRanOrder: Doublep0;
        lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM,lOutImgSum,lOutImgBM,lOutImgT,lOutImgAUC: singlep;
        lSumImg: Smallintp;
        lPlankImg: byteP;
        lVoxPerPlank,lnPlanks,lTotalMemory,lnVoxTested,lVolVox: int64;
        lUniqueOrders,lThread,lThreadStart,lThreadInc,lThreadEnd,
        lPos2,lPosPct,lPos,lPlankImgPos,lPlank,lStartVox,lEndVox: integer;
        lOverlapRA: Overlapp;
     {$IFNDEF FPC} lStartTime :DWord;{$ENDIF}
begin
     {$IFNDEF FPC} lStartTime := GetTickCount;{$ENDIF}
     result := false;
     lSumImg := nil;
     lPlankImg := nil;
     lOutImgSum := nil;
     lOutImgBM := nil;
     lOutImgT := nil;
     lOutImgAUC := nil;
     lOverlapRA := nil;
     lUniqueOrders := 0;
     if lPrefs.Ltest then begin
        lPrefs.Ttest := false;
        lPrefs.BMtest := false;
     end else if (not lPrefs.Ttest) and (not lPrefs.BMtest) then begin//not binomial
        MainForm.NPMmsg('Error no tests specified');
        exit;
     end;
        MainForm.NPMmsg('Permutations = ' +IntToStr(lPrefs.nPermute));
	MainForm.NPMmsg('Analysis began = ' +TimeToStr(Now));
	lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];
	if (lVolVox < 1) then goto 667;
  if not MakeSum( lImages, lMaskHdr, lSumImg) then goto 667;
  lnVoxTested := ThreshSumImg(lSumImg,lVolVox,lPrefs.nCrit);
	MainForm.NPMmsg('Voxels damaged in at least '+inttostr(lPrefs.nCrit)+' individuals = ' +Floattostr(lnVoxTested));
  if lnVoxTested < 1 then begin
	   MainForm.NPMmsg('Error: no voxels damaged in at least '+inttostr(lPrefs.nCrit)+' individuals.');
     goto 667;
  end;
  if (lPrefs.ExplicitMaskName <> '') then begin
    lnVoxTested := ExplicitMaskSumImg (lPrefs.ExplicitMaskName, lSumImg, lVolVox);
	  MainForm.NPMmsg('Voxels also non-zero in mask '+lPrefs.ExplicitMaskName+' = ' +Floattostr(lnVoxTested));
    if lnVoxTested < 1 then begin
	    MainForm.NPMmsg('Error: no remaing voxels also non-zero in mask '+lPrefs.ExplicitMaskName);
      goto 667;
    end;
  end;

        //compute planks and acquire memory
	lTotalMemory := lnVoxTested * lImages.Count;
	if (lTotalMemory = 0)  then goto 667; //no data
	lnPlanks := trunc(lTotalMemory/kPlankSz ) + 1;
	MainForm.NPMmsg('Memory planks = ' +Floattostr(lTotalMemory/kPlankSz));
        if (lnPlanks = 1) then begin
            lVoxPerPlank := lnVoxTested; //we can do this in a single pass
            getmem(lPlankImg,lTotalMemory)
        end else begin
	    getmem(lPlankImg,kPlankSz);
            lVoxPerPlank :=  kPlankSz div lImages.Count;
        end;
        //spatial maps for results
        getmem(lOutImgSum,lVolVox*sizeof(single));
        getmem(lOutImgBM,lVolVox*sizeof(single));
        getmem(lOutImgT,lVolVox*sizeof(single));
        getmem(lOutImgAUC,lVolVox*sizeof(single));
        //initialize memory
        MainForm.InitPermute (lImages.Count, lPrefs.nPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);
	for lPos := 1 to lVolVox do begin
                lOutImgSum^[lPos] := 0;
		lOutImgBM^[lPos] := 0;
		lOutImgT^[lPos] := 0;
		lOutImgAUC^[lPos] := 0;
	end;
        //next create permuted BM bounds
        if lPrefs.BMtest then begin
           MainForm.NPMmsg('Generating BM permutation thresholds');
           MainForm.Refresh;
	   createArray64(lObsp,lObs,lImages.Count);
           for lPos := 1 to lImages.Count do
               lObs^[lPos-1] := lSymptomRA^[lPos];
           genBMsim (lImages.Count, lObs);
           freemem(lObsp);
        end;
        if lPrefs.NULP then
           getmem(lOverlapRA,lnVoxTested* sizeof(TLesionPattern));
        if lPrefs.Ltest then
           ClearThreadDataPvals(gnCPUThreads,lPrefs.nPermute)
        else
            ClearThreadData(gnCPUThreads,lPrefs.nPermute) ;
        //load and process data
	lStartVox := 1;
	lEndVox := 0;
	for lPlank := 1 to lnPlanks do begin
		MainForm.NPMmsg('Computing plank = ' +Inttostr(lPlank)+' of '+inttostr(lnPlanks));
		lEndVox := lEndVox + lVoxPerPlank;
		if lEndVox > lnVoxTested then begin
			lVoxPerPlank := lnVoxTested-lStartVox+1{lVoxPerPlank - (lEndVox-lVolVox)};
			lEndVox := lnVoxTested;
		end;
		lPlankImgPos := 1;
		for lPos := 1 to lImages.Count do begin
			if not LoadImg8Masked(lImages[lPos-1], lPlankImg,lSumImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then
				goto 667;
			lPlankImgPos := lPlankImgPos + lVoxPerPlank;
		end;//for each image
                lThreadStart := 1;
                lThreadInc := lVoxPerPlank  div gnCPUThreads;
                lThreadEnd := lThreadInc;
                Application.processmessages;
                for lThread := 1 to gnCPUThreads do begin
                    if lThread = gnCPUThreads then
                       lThreadEnd := lVoxPerPlank; //avoid integer rounding error
                    if lPrefs.Ltest then begin
                       with TLesionBinom.Create (MainForm.ProgressBar1,false,true,lPrefs.nCrit, lPrefs.nPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,0,lPlankImg,lOutImgSum,lOutImgBM,lOutIMgT{not used},lOutImgAUC,lSymptomRA) do
                         {$IFDEF FPC} OnTerminate := @MainForm.ThreadDone; {$ELSE}OnTerminate := MainForm.ThreadDone;{$ENDIF}
                    end else begin
                        with TLesionContinuous.Create (MainForm.ProgressBar1,lPrefs.ttest,lPrefs.BMtest,lPrefs.nCrit, lPrefs.nPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,0,lPlankImg,lOutImgSum,lOutImgBM,lOutImgT,lOutImgAUC,lSymptomRA) do
                        //with TLesionContinuous.Create (MainForm.ProgressBar1,lttest,lBM,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,lPlankImg,lOutImgSum,lOutImgBM,lOutImgT,lSymptomRA) do
                         {$IFDEF FPC} OnTerminate := @MainForm.ThreadDone; {$ELSE}OnTerminate := MainForm.ThreadDone;{$ENDIF}
                    end;
                    inc(gThreadsRunning);
                    lThreadStart := lThreadEnd + 1;
                    lThreadEnd :=lThreadEnd + lThreadInc;
                end; //for each thread
                repeat
                      Application.processmessages;
                until gThreadsRunning = 0;
                Application.processmessages;
                //end of threading

                if lPrefs.NULP then
                   NULPcount (lPlankImg, lVoxPerPlank,lImages.Count, lUniqueOrders, lOverlapRA);

		lStartVox := lEndVox + 1;
	end;
        //calculate max per thread
        SumThreadData(gnCPUThreads,lPrefs.nPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);

        //data in maps is stored in voxels 1..lnVoxTested - put in spatial order
        reformat(lOutImgSum,lSumImg,lVolVox);
        reformat(lOutImgBM,lSumImg,lVolVox);
        reformat(lOutImgT,lSumImg,lVolVox);
        reformat(lOutImgAUC,lSumImg,lVolVox);
        lThreshBonf := MainForm.reportBonferroni('Std',lnVoxTested);
        if lPrefs.NULP then
           lThreshBonf := MainForm.reportBonferroni('Number of Unique Lesion Patterns',lUniqueOrders);


        //next: save data
	MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);
        //save sum map
        lOutNameMod := ChangeFilePostfixExt(lOutName,'Sum'+lFactName,'.hdr');
        if lPrefs.Run < 1 then
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);
        //save Area Under Curve
        lOutNameMod := ChangeFilePostfixExt(lOutName,'rocAUC'+lFactName,'.hdr');
        if lPrefs.Run < 1 then
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgAUC,1);
        //create new header - subsequent images will use Z-scores
        MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );
        if (lPrefs.Run < 1) and (Sum2Power(lOutImgSum,lVolVox,lImages.Count,lPrefs.nCrit, lPrefs.LTest)) then begin
           lOutNameMod := ChangeFilePostfixExt(lOutName,'Power'+lFactName,'.hdr');
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);
        end;
        //if lPrefs.Run > 0 then //terrible place to do this - RAM problems, but need value to threshold maps
        //   lThreshNULP := MainForm.reportBonferroni('Unique overlap',CountOverlap2 (lImages, lPrefs.nCrit,lnVoxTested,lPlankImg));
        if lPrefs.ttest then begin //save Ttest
           //next: convert t-scores to z scores
           for lPos := 1 to lVolVox do
               lOutImgT^[lPos] := TtoZ (lOutImgT^[lPos],lImages.Count-2);
           for lPos := 1 to lPrefs.nPermute do begin
               lPermuteMaxT^[lPos] := TtoZ (lPermuteMaxT^[lPos],lImages.Count-2);
               lPermuteMinT^[lPos] := TtoZ (lPermuteMinT^[lPos],lImages.Count-2);
           end;
           lThreshFDR := MainForm.reportFDR ('ttest', lVolVox, lnVoxTested, lOutImgT);
           lThreshPermute := MainForm.reportPermute('ttest',lPrefs.nPermute,lPermuteMaxT, lPermuteMinT);
	   lOutNameMod := ChangeFilePostfixExt(lOutName,'ttest'+lFactName,'.hdr');
           {$IFNDEF FPC}
           if lPrefs.Run > 0 then begin
              MainForm.NPMmsgAppend('threshtt,'+inttostr(lPrefs.Run)+','+inttostr(MainForm.ThreshMap(lThreshBonf,lVolVox,lOutImgT))+','+realtostr(lThreshNULP,3)+','+realtostr(lThreshPermute,3)+','+realtostr(lThreshBonf,3)+','+inttostr(round((GetTickCount - lStartTime)/1000)));
           end;
           {$ENDIF}
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgT,1);
        end;
        if lPrefs.LTest then begin
           PtoZpermute (lPrefs.nPermute, lPermuteMaxT, lPermuteMinT);
           lOutNameMod := ChangeFilePostfixExt(lOutName,'L'+lFactName,'.hdr');
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgBM,1);
           MainForm.reportFDR ('L', lVolVox, lnVoxTested, lOutImgBM);
           MainForm.reportPermute('L',lPrefs.nPermute,lPermuteMaxT, lPermuteMinT);
        end;//Liebermeister
        if lPrefs.BMtest then begin //save Brunner Munzel
           lThreshFDR :=  MainForm.reportFDR ('BM', lVolVox, lnVoxTested, lOutImgBM);
           lThreshPermute := MainForm.reportPermute('BM',lPrefs.nPermute,lPermuteMaxBM, lPermuteMinBM);
           lOutNameMod := ChangeFilePostfixExt(lOutName,'BM'+lFactName,'.hdr');
           if lPrefs.Run > 0 then
              MainForm.NPMmsgAppend('threshbm,'+inttostr(lPrefs.Run)+','+inttostr(MainForm.ThreshMap(lThreshBonf,lVolVox,lOutImgBM))+','+realtostr(lThreshNULP,3)+','+realtostr(lThreshPermute,3)+','+realtostr(lThreshBonf,3));
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgBM,1);
        end;
	MainForm.NPMmsg('Analysis finished = ' +TimeToStr(Now));
     {$IFNDEF FPC} MainForm.NPMmsg('Processing Time = ' +inttostr(round((GetTickCount - lStartTime)/1000)));{$ENDIF}

        lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes'+lFactName,'.txt');
        MainForm.MsgSave(lOutNameMod);
        //all done
        result := true;//all done without aborting
667: // free memory and report error
        if lPlankImg <> nil then freemem(lPlankImg);
        if lSumImg <> nil then freemem(lSumImg);
        if lOutImgSum <> nil then freemem(lOutImgSum);
        if lOutImgBM <> nil then freemem(lOutImgBM);
        if lOutImgT <> nil then freemem(lOutImgT);
        if lOutImgAUC <> nil then freemem(lOutImgAUC);
        if lOverlapRA <> nil then freemem(lOverlapRA);

        if not result then
	   MainForm.NPMmsg('Unable to complete analysis.');
        MainForm.ProgressBar1.Position := 0;
end; //TurboLDM

end.