← Back to branch summary

~ubuntu-branches/ubuntu/trusty/mricron/trusty-proposed

« back to all changes in this revision

Viewing changes to npm/old/lesion.pas

  • Committer: Bazaar Package Importer
  • Author(s): Michael Hanke
  • Date: 2010-07-29 22:07:43 UTC
  • Revision ID: james.westby@ubuntu.com-20100729220743-q621ts2zj806gu0n
Tags: upstream-0.20100725.1~dfsg.1
Import upstream version 0.20100725.1~dfsg.1

Show diffs side-by-side

added added

removed removed

Lines of Context:
npm/old/lesion.pas
 
1
unit lesion;
 
2
interface
 
3
{$H+}
 
4
uses
 
5
  define_types,SysUtils,
 
6
part,StatThds,statcr,StatThdsUtil,Brunner,DISTR,nifti_img, hdr,
 
7
   Messages,  Classes, Graphics, Controls, Forms, Dialogs,
 
8
StdCtrls,  ComCtrls,ExtCtrls,Menus, overlap,ReadInt,lesion_pattern,stats,LesionStatThds,nifti_hdr,
 
9
 
 
10
{$IFDEF FPC} LResources,gzio2,
 
11
{$ELSE} gziod,associate,{$ENDIF}   //must be in search path, e.g. C:\pas\mricron\npm\math
 
12
{$IFNDEF UNIX} Windows, {$ENDIF}
 
13
upower,firthThds,firth,IniFiles,cpucount,userdir,math,
 
14
regmult,utypes;
 
15
 
 
16
 
 
17
function LesionNPMAnalyze2 (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lnCrit,lRun,lnPermute: integer; var lSymptomRA: SingleP;var lFactname,lOutName: string; lttest,lBM: boolean): boolean;
 
18
function LesionNPMAnalyzeBinomial2 (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lnCrit,lnPermute: integer; var lSymptomRA: SingleP; var lFactname,lOutName: string): boolean;
 
19
var
 
20
   gNULP,gROI: boolean;
 
21
implementation
 
22
 
 
23
uses npmform;
 
24
 
 
25
{$DEFINE NOTmedianfx}
 
26
function LesionNPMAnalyze2 (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lnCrit,lRun,lnPermute: integer; var lSymptomRA: SingleP;var lFactname,lOutName: string; lttest,lBM: boolean): boolean;
 
27
label
 
28
        123,667;
 
29
var
 
30
        lOutNameMod: string;
 
31
        lPlankImg: byteP;
 
32
        lOutImgSum,lOutImgBM,lOutImgT,lOutImgAUC,
 
33
        lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP;
 
34
        lPos,lPlank,lThread: integer;
 
35
        lVolVox,lMinMask,lMaxMask,lTotalMemory,lnPlanks,lVoxPerPlank,
 
36
        lThreadStart,lThreadEnd,lThreadInc,lnLesion,//,lnPermute,
 
37
        lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lPosPct: int64;
 
38
        lT,lBMz,  lSum,lThresh,lThreshPermute,lThreshBonf,lThreshNULP :double;
 
39
        lObsp: pointer;
 
40
        lObs: Doublep0;
 
41
        lStatHdr: TNIfTIhdr;
 
42
        lFdata: file;
 
43
        lRanOrderp: pointer;
 
44
        lRanOrder: Doublep0;
 
45
        lPlankAllocated: boolean;
 
46
        //lttest,lBM: boolean;
 
47
        {$IFDEF medianfx}
 
48
        lmedianFX,lmeanFX,lsummean,lsummedian: double;
 
49
        lmediancount: integer;
 
50
        {$ENDIF}
 
51
begin
 
52
        //lttest:= ttestmenu.checked;
 
53
        //lBM := BMmenu.checked;
 
54
        lPlankAllocated := false;
 
55
        //lnPermute := MainForm.ReadPermute;
 
56
        MainForm.NPMmsg('Permutations = ' +IntToStr(lnPermute));
 
57
        MainForm.NPMmsg('Analysis began = ' +TimeToStr(Now));
 
58
        lTotalMemory := 0;
 
59
        lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];
 
60
        if (lVolVox < 1) then goto 667;
 
61
        lMinMask := 1;
 
62
        lMaxMask := lVolVox;
 
63
        lVoxPerPlank :=  kPlankSz div lImages.Count div sizeof(byte) ;
 
64
        if (lVoxPerPlank = 0) then goto 667; //no data
 
65
        lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;
 
66
        if (lTotalMemory = 0)  then goto 667; //no data
 
67
        lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;
 
68
        MainForm.NPMmsg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));
 
69
        MainForm.NPMmsg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));
 
70
        if (lnPlanks = 1) then
 
71
            getmem(lPlankImg,lTotalMemory) //assumes 1bpp
 
72
        else
 
73
            getmem(lPlankImg,kPlankSz);
 
74
        lPlankAllocated := true;
 
75
        lStartVox := lMinMask;
 
76
        lEndVox := lMinMask-1;
 
77
        {$IFDEF medianfx}
 
78
        lsummean := 0;
 
79
        lsummedian:= 0;
 
80
        lmediancount := 0;
 
81
        {$ENDIF}
 
82
        for lPos := 1 to lImages.Count do
 
83
                if gScaleRA[lPos] = 0 then
 
84
                        gScaleRA[lPos] := 1;
 
85
        createArray64(lObsp,lObs,lImages.Count);
 
86
        getmem(lOutImgSum,lVolVox* sizeof(single));
 
87
        getmem(lOutImgBM,lVolVox* sizeof(single));
 
88
        getmem(lOutImgT,lVolVox* sizeof(single));
 
89
        getmem(lOutImgAUC,lVolVox* sizeof(single));
 
90
        MainForm.InitPermute (lImages.Count, lnPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);
 
91
        for lPos := 1 to lVolVox do begin
 
92
                lOutImgSum^[lPos] := 0;
 
93
                lOutImgBM^[lPos] := 0;
 
94
                lOutImgT^[lPos] := 0;
 
95
                lOutImgAUC^[lPos] := 0;
 
96
        end;
 
97
        //next create permuted BM bounds
 
98
        if lBM then begin
 
99
           MainForm.NPMmsg('Generating BM permutation thresholds');
 
100
           MainForm.Refresh;
 
101
           for lPos := 1 to lImages.Count do
 
102
               lObs^[lPos-1] := lSymptomRA^[lPos];
 
103
           genBMsim (lImages.Count, lObs);
 
104
        end;
 
105
         ClearThreadData(gnCPUThreads,lnPermute) ;
 
106
        for lPlank := 1 to lnPlanks do begin
 
107
                MainForm.NPMmsg('Computing plank = ' +Inttostr(lPlank));
 
108
        MainForm.Refresh;
 
109
        Application.processmessages;
 
110
                lEndVox := lEndVox + lVoxPerPlank;
 
111
                if lEndVox > lMaxMask then begin
 
112
                        lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);
 
113
                        lEndVox := lMaxMask;
 
114
                end;
 
115
                lPlankImgPos := 1;
 
116
                for lPos := 1 to lImages.Count do begin
 
117
                        if not LoadImg8(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then
 
118
                                goto 667;
 
119
                        lPlankImgPos := lPlankImgPos + lVoxPerPlank;
 
120
                end;//for each image
 
121
                //threading start
 
122
                lThreadStart := 1;
 
123
                lThreadInc := lVoxPerPlank  div gnCPUThreads;
 
124
                lThreadEnd := lThreadInc;
 
125
                Application.processmessages;
 
126
                for lThread := 1 to gnCPUThreads do begin
 
127
                    if lThread = gnCPUThreads then
 
128
                       lThreadEnd := lVoxPerPlank; //avoid integer rounding error
 
129
                    with TLesionContinuous.Create (MainForm.ProgressBar1,lttest,lBM,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,0,lPlankImg,lOutImgSum,lOutImgBM,lOutImgT,lOutImgAUC,lSymptomRA) do
 
130
                    //with TLesionContinuous.Create (MainForm.ProgressBar1,lttest,lBM,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,lPlankImg,lOutImgSum,lOutImgBM,lOutImgT,lSymptomRA) do
 
131
                         {$IFDEF FPC} OnTerminate := @MainForm.ThreadDone; {$ELSE}OnTerminate := MainForm.ThreadDone;{$ENDIF}
 
132
                    inc(gThreadsRunning);
 
133
                    lThreadStart := lThreadEnd + 1;
 
134
                    lThreadEnd :=lThreadEnd + lThreadInc;
 
135
                end; //for each thread
 
136
                repeat
 
137
                      Application.processmessages;
 
138
                until gThreadsRunning = 0;
 
139
                Application.processmessages;
 
140
                //threading end
 
141
                lStartVox := lEndVox + 1;
 
142
        end;
 
143
        //freemem(lPlankImg);
 
144
        //lPlankAllocated := false;
 
145
        lThreshPermute := 0;
 
146
        lnVoxTested := SumThreadData(gnCPUThreads,lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);
 
147
        //next report findings
 
148
        if lnVoxTested < 1 then begin
 
149
           MainForm.NPMmsg('**Error: no voxels tested: no regions lesioned in at least '+inttostr(lnCrit)+' patients**');
 
150
           goto 123;
 
151
        end;
 
152
        MainForm.NPMmsg('Voxels tested = ' +Inttostr(lnVoxTested));
 
153
        {$IFDEF medianfx}
 
154
        MainForm.NPMmsg('Average MEAN effect size = ' +realtostr((lsummean/lmediancount),3));
 
155
        MainForm.NPMmsg('Average MEDIAN effect size = ' +realtostr((lsummedian/lmediancount),3));
 
156
        {$ENDIF}
 
157
        MainForm.NPMmsg('Only tested voxels with more than '+inttostr(lnCrit)+' lesions');
 
158
        //Next: save results from permutation thresholding....
 
159
        lThreshBonf := MainForm.reportBonferroni('Std',lnVoxTested);
 
160
 
 
161
        //next: save data
 
162
        MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);
 
163
//save sum map
 
164
        lOutNameMod := ChangeFilePostfixExt(lOutName,'Sum'+lFactName,'.hdr');
 
165
        if lRun < 1 then
 
166
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);
 
167
        //save Area Under Curve
 
168
        lOutNameMod := ChangeFilePostfixExt(lOutName,'rocAUC'+lFactName,'.hdr');
 
169
        if lRun < 1 then
 
170
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgAUC,1);
 
171
//create new header - subsequent images will use Z-scores
 
172
        MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );
 
173
        if (lRun < 1) and (Sum2PowerCont(lOutImgSum,lVolVox,lImages.Count)) then begin
 
174
           lOutNameMod := ChangeFilePostfixExt(lOutName,'Power'+lFactName,'.hdr');
 
175
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);
 
176
        end;
 
177
        if lRun > 0 then //terrible place to do this - RAM problems, but need value to threshold maps
 
178
           lThreshNULP := MainForm.reportBonferroni('Unique overlap',CountOverlap2 (lImages, lnCrit,lnVoxTested,lPlankImg));
 
179
 
 
180
        //MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );
 
181
if lttest then begin //save Ttest
 
182
        //next: convert t-scores to z scores
 
183
        for lPos := 1 to lVolVox do
 
184
            lOutImgT^[lPos] := TtoZ (lOutImgT^[lPos],lImages.Count-2);
 
185
        for lPos := 1 to lnPermute do begin
 
186
            lPermuteMaxT^[lPos] := TtoZ (lPermuteMaxT^[lPos],lImages.Count-2);
 
187
            lPermuteMinT^[lPos] := TtoZ (lPermuteMinT^[lPos],lImages.Count-2);
 
188
        end;
 
189
        lThresh := MainForm.reportFDR ('ttest', lVolVox, lnVoxTested, lOutImgT);
 
190
        lThreshPermute := MainForm.reportPermute('ttest',lnPermute,lPermuteMaxT, lPermuteMinT);
 
191
        lOutNameMod := ChangeFilePostfixExt(lOutName,'ttest'+lFactName,'.hdr');
 
192
        if lRun > 0 then
 
193
           MainForm.NPMmsgAppend('threshtt,'+inttostr(lRun)+','+inttostr(MainForm.ThreshMap(lThreshNULP,lVolVox,lOutImgT))+','+realtostr(lThreshNULP,3)+','+realtostr(lThreshPermute,3)+','+realtostr(lThreshBonf,3));
 
194
        NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgT,1);
 
195
 
 
196
end;
 
197
if lBM then begin //save Brunner Munzel
 
198
        lThresh :=  MainForm.reportFDR ('BM', lVolVox, lnVoxTested, lOutImgBM);
 
199
        lThreshPermute := MainForm.reportPermute('BM',lnPermute,lPermuteMaxBM, lPermuteMinBM);
 
200
        lOutNameMod := ChangeFilePostfixExt(lOutName,'BM'+lFactName,'.hdr');
 
201
        if lRun > 0 then
 
202
           MainForm.NPMmsgAppend('threshbm,'+inttostr(lRun)+','+inttostr(MainForm.ThreshMap(lThreshNULP,lVolVox,lOutImgBM))+','+realtostr(lThreshNULP,3)+','+realtostr(lThreshPermute,3)+','+realtostr(lThreshBonf,3));
 
203
        NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgBM,1);
 
204
end;
 
205
//next: free dynamic memory
 
206
123:
 
207
        MainForm.FreePermute (lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM,  lRanOrderp);
 
208
        freemem(lOutImgT);
 
209
        freemem(lOutImgAUC);
 
210
        freemem(lOutImgBM);
 
211
        freemem(lOutImgSum);
 
212
        freemem(lObsp);
 
213
        if lPlankAllocated then
 
214
           freemem(lPlankImg);
 
215
        //Next: NULPS - do this after closing all memory - this is a memory hog
 
216
        if gNULP then
 
217
           lThreshNULP := MainForm.reportBonferroni('Unique overlap',CountOverlap (lImages, lnCrit,lnVoxTested));
 
218
        MainForm.NPMmsg('Analysis finished = ' +TimeToStr(Now));
 
219
        lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes'+lFactName,'.txt');
 
220
        MainForm.MsgSave(lOutNameMod);
 
221
        MainForm.ProgressBar1.Position := 0;
 
222
        //if lRun > 0 then
 
223
        //   AX(freeram,freeram,freeram,freeram,freeram,freeram);
 
224
        exit;
 
225
667: //you only get here if you aborted ... free memory and report error
 
226
        if lTotalMemory > 1 then freemem(lPlankImg);
 
227
        MainForm.NPMmsg('Unable to complete analysis.');
 
228
        MainForm.ProgressBar1.Position := 0;
 
229
end; //LesionNPMAnalyze
 
230
 
 
231
function LesionNPMAnalyzeBinomial2 (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lnCrit,lnPermute: integer; var lSymptomRA: SingleP; var lFactname,lOutName: string): boolean;
 
232
label
 
233
     123,667;
 
234
var
 
235
   lVal: single;
 
236
        lOutNameMod: string;
 
237
        lPlankImg: byteP;
 
238
        lOutImgSum,lOutImgL,lOutImgAUC,lDummyImg,
 
239
        lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM: singleP;
 
240
        lPos,lPlank,lThread,lnDeficit: integer;
 
241
        lTotalMemory,lVolVox,lMinMask,lMaxMask,lnPlanks,lVoxPerPlank,
 
242
        lThreadStart,lThreadInc,lThreadEnd, lnLesion,
 
243
        lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lPosPct: int64;
 
244
        lT,  lSum: double;
 
245
        lObsp: pointer;
 
246
        lObs: Doublep0;
 
247
        lStatHdr: TNIfTIhdr;
 
248
        lFdata: file;
 
249
        lRanOrderp: pointer;
 
250
        lRanOrder: Doublep0;
 
251
begin
 
252
        MainForm.NPMmsg('Permutations = ' +IntToStr(lnPermute));
 
253
        //lOutName := lMaskHdr.ImgFileName;
 
254
        //if not SaveHdrName ('Statistical Map', lOutName) then exit;
 
255
        MainForm.NPMmsg('Analysis began = ' +TimeToStr(Now));
 
256
        lTotalMemory := 0;
 
257
        lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];
 
258
        if (lVolVox < 1) then goto 667;
 
259
        lMinMask := 1;
 
260
        lMaxMask := lVolVox;
 
261
        lVoxPerPlank :=  kPlankSz div lImages.Count div sizeof(byte) ;
 
262
        if (lVoxPerPlank = 0) then goto 667; //no data
 
263
        lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;
 
264
        if (lTotalMemory = 0)  then goto 667; //no data
 
265
        lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;
 
266
        MainForm.NPMmsg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));
 
267
        MainForm.NPMmsg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));
 
268
    if (lnPlanks = 1) then
 
269
            getmem(lPlankImg,lTotalMemory) //assumes 1bp
 
270
    else
 
271
            getmem(lPlankImg,kPlankSz);
 
272
        lStartVox := lMinMask;
 
273
        lEndVox := lMinMask-1;
 
274
        for lPos := 1 to lImages.Count do
 
275
                if gScaleRA[lPos] = 0 then
 
276
                        gScaleRA[lPos] := 1;
 
277
        createArray64(lObsp,lObs,lImages.Count);
 
278
        getmem(lOutImgSum,lVolVox* sizeof(single));
 
279
        getmem(lOutImgL,lVolVox* sizeof(single));
 
280
        getmem(lOutImgAUC,lVolVox* sizeof(single));
 
281
 
 
282
        MainForm.InitPermute (lImages.Count, lnPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);
 
283
        for lPos := 1 to lVolVox do begin
 
284
                lOutImgSum^[lPos] := 0;
 
285
                lOutImgL^[lPos] := 0;
 
286
                lOutImgAUC^[lPos] := 0;
 
287
        end;
 
288
        ClearThreadDataPvals(gnCPUThreads,lnPermute) ;
 
289
        for lPlank := 1 to lnPlanks do begin
 
290
            MainForm.ProgressBar1.Position := 1;
 
291
                MainForm.NPMmsg('Computing plank = ' +Inttostr(lPlank));
 
292
                MainForm.Refresh;
 
293
                Application.processmessages;
 
294
                lEndVox := lEndVox + lVoxPerPlank;
 
295
                if lEndVox > lMaxMask then begin
 
296
                        lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);
 
297
                        lEndVox := lMaxMask;
 
298
                end;
 
299
                lPlankImgPos := 1;
 
300
                for lPos := 1 to lImages.Count do begin
 
301
                        if not LoadImg8(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then
 
302
                                goto 667;
 
303
                        lPlankImgPos := lPlankImgPos + lVoxPerPlank;
 
304
                end;//for each image
 
305
                  //threading start
 
306
                lThreadStart := 1;
 
307
                lThreadInc := lVoxPerPlank  div gnCPUThreads;
 
308
                lThreadEnd := lThreadInc;
 
309
                Application.processmessages;
 
310
                for lThread := 1 to gnCPUThreads do begin
 
311
                    if lThread = gnCPUThreads then
 
312
                       lThreadEnd := lVoxPerPlank; //avoid integer rounding error
 
313
                    //with TLesionBinomial.Create (ProgressBar1,false,true,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,666, lDummyImg,lPlankImg,lOutImgSum,lOutImgL,lDummyImg,lSymptomRA) do
 
314
                    with TLesionBinom.Create (MainForm.ProgressBar1,false,true,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,0,lPlankImg,lOutImgSum,lOutImgL,lDummyImg,lOutImgAUC,lSymptomRA) do
 
315
                         {$IFDEF FPC} OnTerminate := @MainForm.ThreadDone; {$ELSE}OnTerminate := MainForm.ThreadDone;{$ENDIF}
 
316
                    inc(gThreadsRunning);
 
317
                    MainForm.NPMmsg('Thread ' +Inttostr(gThreadsRunning)+' = '+inttostr(lThreadStart)+'..'+inttostr(lThreadEnd));
 
318
                    lThreadStart := lThreadEnd + 1;
 
319
                    lThreadEnd :=lThreadEnd + lThreadInc;
 
320
                end; //for each thread
 
321
                repeat
 
322
                      Application.processmessages;
 
323
                until gThreadsRunning = 0;
 
324
                Application.processmessages;
 
325
                //threading end
 
326
                lStartVox := lEndVox + 1;
 
327
        end;
 
328
        lnVoxTested := SumThreadData(gnCPUThreads,lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);
 
329
        for lPos := 1 to lnPermute do begin
 
330
            if (lPermuteMinT^[lPos] > 1.1) or (lPermuteMinT^[lPos] < -1.1) then
 
331
               lPermuteMinT^[lPos] := 0.5;
 
332
            if (lPermuteMaxT^[lPos] > 1.1) or (lPermuteMaxT^[lPos] < -1.1) then
 
333
               lPermuteMaxT^[lPos] := 0.5;
 
334
            lVal := lPermuteMaxT^[lPos];
 
335
            lPermuteMaxT^[lPos] := lPermuteMinT^[lPos];
 
336
            lPermuteMinT^[lPos] := lVal;
 
337
            if lPermuteMaxT^[lPos] < 0 then
 
338
                        lPermuteMaxT^[lPos] := -pNormalInv(abs(lPermuteMaxT^[lPos]))
 
339
            else
 
340
                        lPermuteMaxT^[lPos] := pNormalInv(lPermuteMaxT^[lPos]);
 
341
            if lPermuteMinT^[lPos] < 0 then
 
342
                        lPermuteMinT^[lPos] := -pNormalInv(abs(lPermuteMinT^[lPos]))
 
343
            else
 
344
                        lPermuteMinT^[lPos] := pNormalInv(lPermuteMinT^[lPos]);
 
345
        end;
 
346
 
 
347
 
 
348
 
 
349
        if lnVoxTested < 1 then begin
 
350
           MainForm.NPMmsg('**Error: no voxels tested: no regions lesioned in at least '+inttostr(lnCrit)+' patients**');
 
351
           goto 123;
 
352
        end;
 
353
        //next report findings
 
354
        MainForm.NPMmsg('Voxels tested = ' +Inttostr(lnVoxTested));
 
355
        MainForm.NPMmsg('Only tested voxels with more than '+inttostr(lnCrit)+' lesions');
 
356
        //Next: save results from permutation thresholding....
 
357
        MainForm.reportBonferroni('Std',lnVoxTested);
 
358
 
 
359
        //next: save data
 
360
//savedata
 
361
        MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);
 
362
//save sum map
 
363
        lOutNameMod := ChangeFilePostfixExt(lOutName,'Sum'+lFactName,'.hdr');
 
364
        NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);
 
365
//save Area Under Curve
 
366
           lOutNameMod := ChangeFilePostfixExt(lOutName,'rocAUC'+lFactName,'.hdr');
 
367
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgAUC,1);
 
368
 
 
369
//future images will store Z-scores...
 
370
        MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );
 
371
//save power map
 
372
        lnDeficit := 0;
 
373
        for lPos := 1 to lImages.Count do
 
374
            if lSymptomRA^[lPos] = 0 then
 
375
               inc(lnDeficit);
 
376
        if Sum2PowerBinom(lOutImgSum,lVolVox,lImages.Count,lnDeficit) then begin
 
377
           lOutNameMod := ChangeFilePostfixExt(lOutName,'Power'+lFactName,'.hdr');
 
378
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);
 
379
        end;
 
380
        //save Liebermeister
 
381
 
 
382
        lOutNameMod := ChangeFilePostfixExt(lOutName,'L'+lFactName,'.hdr');
 
383
        NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgL,1);
 
384
        //save end
 
385
       MainForm.reportFDR ('L', lVolVox, lnVoxTested, lOutImgL);
 
386
       MainForm.reportPermute('L',lnPermute,lPermuteMaxT, lPermuteMinT);
 
387
 
 
388
123:
 
389
//next: free dynamic memory
 
390
        MainForm.FreePermute (lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM,  lRanOrderp);
 
391
        freemem(lOutImgL);
 
392
        freemem(lOutImgAUC);
 
393
        freemem(lOutImgSum);
 
394
        freemem(lObsp);
 
395
        freemem(lPlankImg);
 
396
        //Next: NULPS  - do this at the end, it is a memory hog!
 
397
        if gNULP then
 
398
           MainForm.reportBonferroni('Unique overlap',CountOverlap (lImages, lnCrit,lnVoxTested));
 
399
 
 
400
        MainForm.NPMmsg('Analysis finished = ' +TimeToStr(Now));
 
401
        lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes'+lFactName,'.txt');
 
402
        MainForm.MsgSave(lOutNameMod);
 
403
 
 
404
        MainForm.ProgressBar1.Position := 0;
 
405
        exit;
 
406
667: //you only get here if you aborted ... free memory and report error
 
407
        if lTotalMemory > 1 then freemem(lPlankImg);
 
408
        MainForm.NPMmsg('Unable to complete analysis.');
 
409
        MainForm.ProgressBar1.Position := 0;
 
410
end;
 
411
 
 
412
 
 
413
 
 
414
 
 
415
 
 
416
end.