← Back to branch summary

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

« back to all changes in this revision

Viewing changes to npm/xanacom.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/xanacom.pas
 
1
unit anacom;
 
2
interface
 
3
{$H+}
 
4
uses
 
5
  define_types,SysUtils,part,StatThds,statcr,StatThdsUtil,Brunner,
 
6
  DISTR,nifti_img, hdr,filename,Messages,  Classes, Graphics,
 
7
  Controls, Forms, Dialogs,StdCtrls,ComCtrls,ExtCtrls,Menus, overlap,
 
8
  ReadInt,lesion_pattern,stats,LesionStatThds,nifti_hdr,
 
9
  upower,firthThds,firth,IniFiles,cpucount,userdir,math,
 
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
  regmult,utypes;
 
14
 
 
15
  function AnacomLesionNPMAnalyze (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lnCrit,lRun,lnControl: integer; var lSymptomRA,lControlSymptomRA: SingleP;var lFactname,lOutName: string; lttestIn,lBMIn: boolean): boolean;
 
16
  procedure DoAnaCOM;
 
17
  function readTxt (lFilename: string; var lnObservations : integer; var ldataRA1: singlep): boolean;
 
18
 
 
19
 
 
20
implementation
 
21
 
 
22
uses npmform;
 
23
 
 
24
{$DEFINE NOTmedianfx}
 
25
function AnacomLesionNPMAnalyze (var lImages: TStrings; var lMaskHdr: TMRIcroHdr; lnCrit,lRun,lnControl: integer; var lSymptomRA,lControlSymptomRA: SingleP;var lFactname,lOutName: string; lttestIn,lBMIn: boolean): boolean;
 
26
label
 
27
        123,667;
 
28
var
 
29
        lOutNameMod: string;
 
30
        lPlankImg: byteP;
 
31
        lOutImgSum,lOutImgBM,lOutImgT,
 
32
        lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM,lCombinedSymptomRA: singleP;
 
33
        lPos,lPlank,lThread,lnControlsPlusPatients: integer;
 
34
        lVolVox,lMinMask,lMaxMask,lTotalMemory,lnPlanks,lVoxPerPlank,
 
35
        lThreadStart,lThreadEnd,lThreadInc,lnLesion,lnPermute,
 
36
        lPos2,lPos2Offset,lStartVox,lEndVox,lPlankImgPos,lnTests,lnVoxTested,lPosPct: int64;
 
37
        lT,lBMz,  lSum,lThresh,lThreshBonf,lThreshPermute,lThreshNULP :double;
 
38
        lObsp: pointer;
 
39
        lObs: Doublep0;
 
40
        lStatHdr: TNIfTIhdr;
 
41
        lFdata: file;
 
42
        lRanOrderp: pointer;
 
43
        lRanOrder: Doublep0;
 
44
        lSave,lBM,lttest,lLtest: boolean;
 
45
        lnControlNeg: integer;
 
46
 
 
47
        {$IFDEF medianfx}
 
48
        lmedianFX,lmeanFX,lsummean,lsummedian: double;
 
49
        lmediancount: integer;
 
50
        {$ENDIF}
 
51
begin
 
52
        lSave := true;
 
53
        lnControlNeg := lnControl; //negative for binomial test
 
54
        lttest := lttestin;
 
55
        lbm := lbmin;
 
56
        if (not (lttest)) and (not (lbm)) then begin
 
57
           lLtest := true;
 
58
           lBM := true;
 
59
           lnControlNeg := -lnControl;
 
60
        end;
 
61
        //lttest:= ttestmenu.checked;
 
62
        //lBM := BMmenu.checked;
 
63
        if lnControl < 1 then begin
 
64
           MainForm.NPMmsg('AnaCOM aborted - need data from at least 1 control individual');
 
65
           exit;
 
66
        end;
 
67
        lnPermute := 0;//MainForm.ReadPermute;
 
68
        MainForm.NPMmsg('Permutations = ' +IntToStr(lnPermute));
 
69
        MainForm.NPMmsg('Analysis began = ' +TimeToStr(Now));
 
70
        lTotalMemory := 0;
 
71
        lVolVox := lMaskHdr.NIFTIhdr.dim[1]*lMaskHdr.NIFTIhdr.dim[2]* lMaskHdr.NIFTIhdr.dim[3];
 
72
        if (lVolVox < 1) then goto 667;
 
73
        lMinMask := 1;
 
74
        lMaxMask := lVolVox;
 
75
        lVoxPerPlank :=  kPlankSz div lImages.Count div sizeof(byte) ;
 
76
        if (lVoxPerPlank = 0) then goto 667; //no data
 
77
        lTotalMemory := ((lMaxMask+1)-lMinMask) * lImages.Count;
 
78
        if (lTotalMemory = 0)  then goto 667; //no data
 
79
        lnPlanks := trunc(lTotalMemory/(lVoxPerPlank*lImages.Count) ) + 1;
 
80
        MainForm.NPMmsg('Memory planks = ' +Floattostr(lTotalMemory/(lVoxPerPlank*lImages.Count)));
 
81
        MainForm.NPMmsg('Max voxels per Plank = ' +Floattostr(lVoxPerPlank));
 
82
        if (lnPlanks = 1) then
 
83
            getmem(lPlankImg,lTotalMemory) //assumes 1bpp
 
84
        else
 
85
            getmem(lPlankImg,kPlankSz);
 
86
        lStartVox := lMinMask;
 
87
        lEndVox := lMinMask-1;
 
88
        {$IFDEF medianfx}
 
89
        lsummean := 0;
 
90
        lsummedian:= 0;
 
91
        lmediancount := 0;
 
92
        {$ENDIF}
 
93
        for lPos := 1 to lImages.Count do
 
94
                if gScaleRA[lPos] = 0 then
 
95
                        gScaleRA[lPos] := 1;
 
96
        lnControlsPlusPatients := lImages.Count+lnControl;
 
97
        createArray64(lObsp,lObs,lnControlsPlusPatients);
 
98
        getmem(lOutImgSum,lVolVox* sizeof(single));
 
99
        getmem(lOutImgBM,lVolVox* sizeof(single));
 
100
        getmem(lOutImgT,lVolVox* sizeof(single));
 
101
        MainForm.InitPermute (lImages.Count, lnPermute, lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM, lRanOrderp, lRanOrder);
 
102
        for lPos := 1 to lVolVox do begin
 
103
                lOutImgSum^[lPos] := 0;
 
104
                lOutImgBM^[lPos] := 0;
 
105
                lOutImgT^[lPos] := 0;
 
106
        end;
 
107
        //sumptom array for lesions AND controls
 
108
        for lPos := 1 to lImages.Count do
 
109
            lObs^[lPos-1] := lSymptomRA^[lPos];
 
110
        for lPos := 1 to lnControl do
 
111
            lObs^[lPos-1+lImages.Count] :=  lControlSymptomRA^[lPos];
 
112
        getmem(lCombinedSymptomRA,lnControlsPlusPatients* sizeof(single));
 
113
        for lPos := 1 to lnControlsPlusPatients do
 
114
            lCombinedSymptomRA^[lPos] := lObs^[lPos-1];
 
115
        //next create permuted BM bounds
 
116
        if lBM then begin
 
117
           MainForm.NPMmsg('Generating BM permutation thresholds');
 
118
           MainForm.Refresh;
 
119
           //for lPos := 1 to lImages.Count do
 
120
           //    lObs^[lPos-1] := lSymptomRA^[lPos];
 
121
           genBMsim (lnControlsPlusPatients, lObs);
 
122
        end;
 
123
         ClearThreadData(gnCPUThreads,lnPermute) ;
 
124
        for lPlank := 1 to lnPlanks do begin
 
125
                MainForm.NPMmsg('Computing plank = ' +Inttostr(lPlank));
 
126
        MainForm.Refresh;
 
127
        Application.processmessages;
 
128
                lEndVox := lEndVox + lVoxPerPlank;
 
129
                if lEndVox > lMaxMask then begin
 
130
                        lVoxPerPlank := lVoxPerPlank - (lEndVox-lMaxMask);
 
131
                        lEndVox := lMaxMask;
 
132
                end;
 
133
                lPlankImgPos := 1;
 
134
                for lPos := 1 to lImages.Count do begin
 
135
                        if not LoadImg8(lImages[lPos-1], lPlankImg, lStartVox, lEndVox,round(gOffsetRA[lPos]),lPlankImgPos,gDataTypeRA[lPos],lVolVox) then
 
136
                                goto 667;
 
137
                        lPlankImgPos := lPlankImgPos + lVoxPerPlank;
 
138
                end;//for each image
 
139
                //threading start
 
140
                lThreadStart := 1;
 
141
                lThreadInc := lVoxPerPlank  div gnCPUThreads;
 
142
                lThreadEnd := lThreadInc;
 
143
                Application.processmessages;
 
144
                for lThread := 1 to gnCPUThreads do begin
 
145
                    if lThread = gnCPUThreads then
 
146
                       lThreadEnd := lVoxPerPlank; //avoid integer rounding error
 
147
 
 
148
                        with TLesionContinuous.Create (MainForm.ProgressBar1,lttest,lBM,lnCrit, lnPermute,lThread,lThreadStart,lThreadEnd,lStartVox,lVoxPerPlank,lImages.Count,lnControlNeg,lPlankImg,lOutImgSum,lOutImgBM,lOutImgT,nil,lCombinedSymptomRA) do
 
149
                         {$IFDEF FPC} OnTerminate := @MainForm.ThreadDone; {$ELSE}OnTerminate := MainForm.ThreadDone;{$ENDIF}
 
150
                    inc(gThreadsRunning);
 
151
                    lThreadStart := lThreadEnd + 1;
 
152
                    lThreadEnd :=lThreadEnd + lThreadInc;
 
153
                end; //for each thread
 
154
                repeat
 
155
                      Application.processmessages;
 
156
                until gThreadsRunning = 0;
 
157
                Application.processmessages;
 
158
                //threading end
 
159
                lStartVox := lEndVox + 1;
 
160
        end;
 
161
                lThreshPermute := 0;
 
162
        lnVoxTested := SumThreadData(gnCPUThreads,lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM);
 
163
        //next report findings
 
164
        if lnVoxTested < 1 then begin
 
165
           MainForm.NPMmsg('**Error: no voxels tested: no regions lesioned in at least '+inttostr(lnCrit)+' patients**');
 
166
           goto 123;
 
167
        end;
 
168
 
 
169
        MainForm.NPMmsg('Voxels tested = ' +Inttostr(lnVoxTested));
 
170
        {$IFDEF medianfx}
 
171
        MainForm.NPMmsg('Average MEAN effect size = ' +realtostr((lsummean/lmediancount),3));
 
172
        MainForm.NPMmsg('Average MEDIAN effect size = ' +realtostr((lsummedian/lmediancount),3));
 
173
        {$ENDIF}
 
174
        MainForm.NPMmsg('Only tested voxels with more than '+inttostr(lnCrit)+' lesions');
 
175
        //Next: save results from permutation thresholding....
 
176
        //Next: save results from permutation thresholding....
 
177
        lThreshBonf := MainForm.reportBonferroni('Std',lnVoxTested);
 
178
        //Next: NULPS
 
179
        if lRun > 0 then //terrible place to do this - RAM problems, but need value to threshold maps
 
180
           lThreshNULP := MainForm.reportBonferroni('Unique overlap',CountOverlap2 (lImages, lnCrit,lnVoxTested,lPlankImg));
 
181
 
 
182
        //lThreshNULP := MainForm.reportBonferroni('Unique overlap',CountOverlap (lImages, lnCrit));
 
183
        //next: save data
 
184
        MakeHdr (lMaskHdr.NIFTIhdr,lStatHdr);
 
185
//save sum map
 
186
        lOutNameMod := ChangeFilePostfixExt(lOutName,'Sum'+lFactName,'.hdr');
 
187
        if (lSave) and (lRun < 1) then
 
188
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);
 
189
//create new header - subsequent images will use Z-scores
 
190
        MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );
 
191
        if (lSave) and (lRun < 1) and (Sum2PowerCont(lOutImgSum,lVolVox,lImages.Count)) then begin
 
192
           lOutNameMod := ChangeFilePostfixExt(lOutName,'Power'+lFactName,'.hdr');
 
193
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgSum,1);
 
194
        end;
 
195
 
 
196
        //MakeStatHdr (lMaskHdr.NIFTIhdr,lStatHdr,-6, 6,1{df},0,lnVoxTested,kNIFTI_INTENT_ZSCORE,inttostr(lnVoxTested) );
 
197
if lttest then begin //save Ttest
 
198
        //next: convert t-scores to z scores
 
199
 
 
200
        if lnControl < 1 then //do not convert t-scores for anaCOM - numbers vary from voxel to voxel...
 
201
           for lPos := 1 to lVolVox do
 
202
               lOutImgT^[lPos] := TtoZ (lOutImgT^[lPos],lImages.Count-2);
 
203
        for lPos := 1 to lnPermute do begin
 
204
            lPermuteMaxT^[lPos] := TtoZ (lPermuteMaxT^[lPos],lImages.Count-2);
 
205
            lPermuteMinT^[lPos] := TtoZ (lPermuteMinT^[lPos],lImages.Count-2);
 
206
        end;
 
207
        lThresh := MainForm.reportFDR ('ttest', lVolVox, lnVoxTested, lOutImgT);
 
208
        lThreshPermute := MainForm.reportPermute('attest',lnPermute,lPermuteMaxT, lPermuteMinT);
 
209
        lOutNameMod := ChangeFilePostfixExt(lOutName,'attest'+lFactName,'.hdr');
 
210
        if lRun > 0 then
 
211
           MainForm.NPMmsgAppend('AnaComthreshtt,'+inttostr(lRun)+','+inttostr(MainForm.ThreshMap(lThreshNULP,lVolVox,lOutImgT))+','+realtostr(lThreshNULP,3)+','+realtostr(lThreshPermute,3)+','+realtostr(lThreshBonf,3));
 
212
        if lSave then
 
213
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgT,1);
 
214
 
 
215
end;
 
216
if lBM then begin //save Mann Whitney
 
217
        lThresh :=  MainForm.reportFDR ('BM', lVolVox, lnVoxTested, lOutImgBM);
 
218
        lThreshPermute := MainForm.reportPermute('aBM',lnPermute,lPermuteMaxBM, lPermuteMinBM);
 
219
        lOutNameMod := ChangeFilePostfixExt(lOutName,'aBM'+lFactName,'.hdr');
 
220
        if lRun > 0 then
 
221
           MainForm.NPMmsgAppend('AnaCOMthreshbm,'+inttostr(lRun)+','+inttostr(MainForm.ThreshMap(lThreshNULP,lVolVox,lOutImgBM))+','+realtostr(lThreshNULP,3)+','+realtostr(lThreshPermute,3)+','+realtostr(lThreshBonf,3));
 
222
        if lSave then
 
223
           NIFTIhdr_SaveHdrImg(lOutNameMod,lStatHdr,true,not IsNifTiMagic(lMaskHdr.NIFTIhdr),true,lOutImgBM,1);
 
224
end;
 
225
//next: free dynamic memory
 
226
123:
 
227
        MainForm.FreePermute (lnPermute,lPermuteMaxT, lPermuteMinT,lPermuteMaxBM, lPermuteMinBM,  lRanOrderp);
 
228
        freemem(lOutImgT);
 
229
        freemem(lOutImgBM);
 
230
        freemem(lOutImgSum);
 
231
        freemem(lObsp);
 
232
        freemem(lPlankImg);
 
233
        freemem(lCombinedSymptomRA);
 
234
        MainForm.NPMmsg('Analysis finished = ' +TimeToStr(Now));
 
235
        lOutNameMod := ChangeFilePostfixExt(lOutName,'Notes'+lFactName,'.txt');
 
236
        if lSave then
 
237
           MainForm.MsgSave(lOutNameMod);
 
238
        MainForm.ProgressBar1.Position := 0;
 
239
        exit;
 
240
667: //you only get here if you aborted ... free memory and report error
 
241
        if lTotalMemory > 1 then freemem(lPlankImg);
 
242
        MainForm.NPMmsg('Unable to complete analysis.');
 
243
        MainForm.ProgressBar1.Position := 0;
 
244
end; //LesionNPMAnalyze
 
245
 
 
246
 
 
247
 
 
248
(*function readCSV2 (lFilename: string; lCol1,lCol2: integer;  var lnObservations : integer; var ldataRA1,ldataRA2: singlep): boolean;
 
249
const
 
250
     kHdrRow = 0;//1;
 
251
     kHdrCol = 0;//1;
 
252
var
 
253
   lNumStr: string;
 
254
   F: TextFile;
 
255
   lTempFloat: double;
 
256
   lCh: char;
 
257
   lnFactors,MaxC,R,C:integer;
 
258
   lError: boolean;
 
259
 
 
260
begin
 
261
     lError := false;
 
262
     result := false;
 
263
         if not fileexists(lFilename) then begin
 
264
            showmessage('Can not find '+lFilename);
 
265
            exit;
 
266
         end;
 
267
         AssignFile(F, lFilename);
 
268
         FileMode := 0;  //Set file access to read only
 
269
         //First pass: determine column height/width
 
270
         Reset(F);
 
271
         C := 0;
 
272
         MaxC := 0;
 
273
         R := 0;
 
274
         while not Eof(F) do begin
 
275
                //read next line
 
276
                //read next line
 
277
                Read(F, lCh);
 
278
                if lCh = '#' then
 
279
                   while not (lCh in [#10,#13]) do
 
280
                           Read(F, lCh)
 
281
                else if not (lCh in [#10,#13,#9,',']) then begin
 
282
                   lNumStr := lNumStr + lCh;
 
283
                end else if lNumStr <> '' then begin
 
284
                        lNumStr := '';
 
285
                        inc(C);
 
286
                        if C > MaxC then
 
287
                                MaxC := C;
 
288
                        if (lCh in [#10,#13]) then begin
 
289
                                C := 0;
 
290
                           inc(R);
 
291
 
 
292
                        end; //eoln
 
293
                end; //if lNumStr <> '' and not tab
 
294
         end;
 
295
         if lNumStr <> '' then  //july06- read data immediately prior to EOF
 
296
        inc(R);
 
297
 
 
298
     if (R <= (kHdrRow+1)) or (MaxC < (kHdrCol+lCol1)) or (MaxC < (kHdrCol+lCol2)) then begin
 
299
         showmessage('problems reading CSV - not enough columns/rows '+inttostr(lCol1)+'  '+inttostr(lCol2));
 
300
         exit;
 
301
     end;
 
302
 
 
303
     lnObservations := R -kHdrRow ; //-1: first row is header....
 
304
     lnFactors := MaxC-1;// -1: first column is Y values
 
305
     //fx(lnObservations,lnFactors);
 
306
 
 
307
     //exit;
 
308
     getmem(ldataRA1,lnObservations*sizeof(single));
 
309
     getmem(ldataRA2,lnObservations*sizeof(single));
 
310
 
 
311
     //second pass
 
312
         Reset(F);
 
313
         C := 1;
 
314
         MaxC := 0;
 
315
         R := 1;
 
316
         lNumStr := '';
 
317
     lTempfloat := 0;
 
318
         while not Eof(F) do begin
 
319
                //read next line
 
320
                Read(F, lCh);
 
321
                if lCh = '#' then
 
322
                   while not (lCh in [#10,#13]) do
 
323
                           Read(F, lCh)
 
324
                else if not (lCh in [#10,#13,#9,',']) then begin
 
325
                   lNumStr := lNumStr + lCh;
 
326
                end else if lNumStr <> '' then begin
 
327
            if (R > kHdrRow) and (C > kHdrCol) then begin
 
328
             if ((C-kHdrCol) = lCol1) or ((C-kHdrCol) = lCol2) then begin
 
329
               if lNumStr = '-' then begin
 
330
                  lTempFloat := 0;
 
331
               end else begin //number
 
332
                try
 
333
                   lTempFloat := strtofloat(lNumStr);
 
334
                except
 
335
                                    on EConvertError do begin
 
336
                                       if not lError then
 
337
                                          showmessage('Empty cells? Error reading CSV file row:'+inttostr(R)+' col:'+inttostr(C)+' - Unable to convert the string '+lNumStr+' to a number');
 
338
                                       lError := true;
 
339
                                       lTempFloat := nan;
 
340
                                    end;
 
341
                end;//except
 
342
                //showmessage(lNumStr);
 
343
                if (C-kHdrCol) = lCol1 then
 
344
                   ldataRA1^[R-kHdrRow] := lTempFloat
 
345
                else if (C-kHdrCol) = lCol2 then
 
346
                    ldataRA2^[R-kHdrRow] := lTempFloat;
 
347
               end; //number
 
348
             end; //col1 or col2
 
349
            end;// else //R > 1
 
350
 
 
351
                        lNumStr := '';
 
352
                        inc(C);
 
353
                        if C > MaxC then
 
354
                                MaxC := C;
 
355
                        if (lCh in [#10,#13]) then begin
 
356
                                C := 1;
 
357
                           inc(R);
 
358
                        end; //eoln
 
359
                end; //if lNumStr <> '' and not tab
 
360
         end;
 
361
     if (lNumStr <> '') and (C = lnFactors) then begin  //unterminated string
 
362
        try
 
363
           lTempFloat := strtofloat(lNumStr);
 
364
        except
 
365
              on EConvertError do begin
 
366
                 if not lError then
 
367
                    showmessage('Empty cells? Error reading CSV file row:'+inttostr(R)+' col:'+inttostr(C)+' - Unable to convert the string '+lNumStr+' to a number');
 
368
                 lError := true;
 
369
                 lTempFloat := nan;
 
370
              end;
 
371
        end;//except
 
372
        ldataRA2^[R-1] := lTempFloat;
 
373
     end;//unterminated string
 
374
     //read finel item
 
375
         CloseFile(F);
 
376
         FileMode := 2;  //Set file access to read/write
 
377
     result := true;
 
378
end;  *)
 
379
 
 
380
function readTxt (lFilename: string; var lnObservations : integer; var ldataRA1: singlep): boolean;
 
381
const
 
382
     kHdrRow = 0;//1;
 
383
     kHdrCol = 0;//1;
 
384
var
 
385
   lCol1: integer;
 
386
   lNumStr: string;
 
387
   F: TextFile;
 
388
   lTempFloat: double;
 
389
   lCh: char;
 
390
   lnFactors,MaxC,R,C:integer;
 
391
   lError: boolean;
 
392
begin
 
393
     lCol1:= 1;
 
394
     lError := false;
 
395
     result := false;
 
396
         if not fileexists(lFilename) then begin
 
397
            showmessage('Can not find '+lFilename);
 
398
            exit;
 
399
         end;
 
400
         AssignFile(F, lFilename);
 
401
         FileMode := 0;  //Set file access to read only
 
402
         //First pass: determine column height/width
 
403
         Reset(F);
 
404
         C := 0;
 
405
         MaxC := 0;
 
406
         R := 0;
 
407
         while not Eof(F) do begin
 
408
                //read next line
 
409
                //read next line
 
410
                Read(F, lCh);
 
411
                if lCh = '#' then
 
412
                   while not (lCh in [#10,#13]) do
 
413
                           Read(F, lCh)
 
414
                else if not (lCh in [#10,#13,#9,',']) then begin
 
415
                   lNumStr := lNumStr + lCh;
 
416
                end else if lNumStr <> '' then begin
 
417
                        lNumStr := '';
 
418
                        inc(C);
 
419
                        if C > MaxC then
 
420
                                MaxC := C;
 
421
                        if (lCh in [#10,#13]) then begin
 
422
                                C := 0;
 
423
                           inc(R);
 
424
 
 
425
                        end; //eoln
 
426
                end; //if lNumStr <> '' and not tab
 
427
         end;
 
428
         if lNumStr <> '' then  //july06- read data immediately prior to EOF
 
429
        inc(R);
 
430
 
 
431
     if (R <= (kHdrRow+1)) or (MaxC < (kHdrCol+lCol1)) then begin
 
432
         showmessage('problems reading CSV - not enough columns/rows ');
 
433
         exit;
 
434
     end;
 
435
 
 
436
     lnObservations := R -kHdrRow ; //-1: first row is header....
 
437
     lnFactors := kHdrCol+lCol1;// -1: first column is Y values
 
438
     //fx(lnObservations,lnFactors);
 
439
 
 
440
     //exit;
 
441
     getmem(ldataRA1,lnObservations*sizeof(single));
 
442
 
 
443
     //second pass
 
444
         Reset(F);
 
445
         C := 1;
 
446
         MaxC := 0;
 
447
         R := 1;
 
448
         lNumStr := '';
 
449
     lTempfloat := 0;
 
450
         while not Eof(F) do begin
 
451
                //read next line
 
452
                Read(F, lCh);
 
453
                if lCh = '#' then
 
454
                   while not (lCh in [#10,#13]) do
 
455
                           Read(F, lCh)
 
456
                else if not (lCh in [#10,#13,#9,',']) then begin
 
457
                   lNumStr := lNumStr + lCh;
 
458
                end else if lNumStr <> '' then begin
 
459
            if (R > kHdrRow) and (C > kHdrCol) then begin
 
460
             if ((C-kHdrCol) = lCol1) {or ((C-kHdrCol) = lCol2)} then begin
 
461
               if lNumStr = '-' then begin
 
462
                  lTempFloat := 0;
 
463
               end else begin //number
 
464
                try
 
465
                   lTempFloat := strtofloat(lNumStr);
 
466
                except
 
467
                                    on EConvertError do begin
 
468
                                       if not lError then
 
469
                                          showmessage('Empty cells? Error reading CSV file row:'+inttostr(R)+' col:'+inttostr(C)+' - Unable to convert the string '+lNumStr+' to a number');
 
470
                                       lError := true;
 
471
                                       lTempFloat := nan;
 
472
                                    end;
 
473
                end;//except
 
474
                //showmessage(lNumStr);
 
475
                if (C-kHdrCol) = lCol1 then begin
 
476
                    //showmessage(lNumStr);
 
477
                   ldataRA1^[R-kHdrRow] := lTempFloat;
 
478
                end;
 
479
                {else if (C-kHdrCol) = lCol2 then
 
480
                    ldataRA2^[R-kHdrRow] := lTempFloat;}
 
481
               end; //number
 
482
             end; //col1 or col2
 
483
            end;// else //R > 1
 
484
 
 
485
                        lNumStr := '';
 
486
                        inc(C);
 
487
                        if C > MaxC then
 
488
                                MaxC := C;
 
489
                        if (lCh in [#10,#13]) then begin
 
490
                                C := 1;
 
491
                           inc(R);
 
492
                        end; //eoln
 
493
                end; //if lNumStr <> '' and not tab
 
494
         end;
 
495
     //showmessage(lNumStr+'  '+inttostr(lnFactors)+'  '+inttostr(C));
 
496
     if (lNumStr <> '') and (C = lnFactors) then begin  //unterminated string
 
497
 
 
498
        try
 
499
           lTempFloat := strtofloat(lNumStr);
 
500
        except
 
501
              on EConvertError do begin
 
502
                 if not lError then
 
503
                    showmessage('Empty cells? Error reading CSV file row:'+inttostr(R)+' col:'+inttostr(C)+' - Unable to convert the string '+lNumStr+' to a number');
 
504
                 lError := true;
 
505
                 lTempFloat := nan;
 
506
              end;
 
507
        end;//except
 
508
         //showmessage(inttostr(R)+'  '+floattostr(lTempFLoat));
 
509
        ldataRA1^[R] := lTempFloat;
 
510
     end;//unterminated string
 
511
     //read finel item
 
512
         CloseFile(F);
 
513
         FileMode := 2;  //Set file access to read/write
 
514
     result := not lError;
 
515
end;
 
516
 
 
517
procedure DoAnaCOM;
 
518
label
 
519
        666;
 
520
var
 
521
   lControlFilename: string;
 
522
   lI, lnControlObservations : integer;
 
523
   lControldata: singlep;
 
524
        //lBinomial: boolean;
 
525
        lFact,lnFactors,lSubj,lnSubj,lnSubjAll,lMaskVoxels,lnCrit: integer;
 
526
        lImageNames,lImageNamesAll:  TStrings;
 
527
        lPredictorList: TStringList;
 
528
        lTemp4D,lMaskname,lOutName,lFactname: string;
 
529
        lMaskHdr: TMRIcroHdr;
 
530
        lMultiSymptomRA,lSymptomRA: singleP;
 
531
begin
 
532
     npmform.MainForm.memo1.lines.clear;
 
533
     npmform.MainForm.memo1.lines.add('AnaCOM analysis requires TXT/CSV format text file.');
 
534
     npmform.MainForm.memo1.lines.add('One row per control participant.');
 
535
     npmform.MainForm.memo1.lines.add('First column is performance of that participant.');
 
536
     npmform.MainForm.memo1.lines.add('Example file:');
 
537
     npmform.MainForm.memo1.lines.add('11');
 
538
     npmform.MainForm.memo1.lines.add('19');
 
539
     npmform.MainForm.memo1.lines.add('2');
 
540
     npmform.MainForm.memo1.lines.add('22');
 
541
     npmform.MainForm.memo1.lines.add('19');
 
542
     npmform.MainForm.memo1.lines.add('6');
 
543
     if not MainForm.OpenDialogExecute('Select text file',false,false,'Text file (*.txt)|*.txt;*.csv') then begin
 
544
           showmessage('AnaCOM aborted: Control data file selection failed.');
 
545
           exit;
 
546
     end; //if not selected
 
547
     lControlFilename := MainForm.OpenHdrDlg.Filename;
 
548
     if (not readTxt (lControlFilename, lnControlObservations,lControldata)) or (lnControlObservations < 1) then begin
 
549
      showmessage('Error reading file '+lControlFilename);
 
550
      exit;
 
551
     end;
 
552
     npmform.MainForm.memo1.lines.add('Control (n='+inttostr(lnControlObservations)+')performance   ['+lControlFilename+']');
 
553
     for lI := 1 to lnControlObservations do
 
554
        npmform.MainForm.memo1.lines.add(inttostr(lI)+' '+floattostr(lControldata^[lI]));
 
555
     //begin - copy
 
556
     lImageNamesAll:= TStringList.Create; //not sure why TStrings.Create does not work???
 
557
     lImageNames:= TStringList.Create; //not sure why TStrings.Create does not work???
 
558
     //next, get 1st group
 
559
     if not MainForm.GetValX(lnSubjAll,lnFactors,lMultiSymptomRA,lImageNamesAll,lnCrit,{,binom}lPredictorList) then begin
 
560
        showmessage('Error with VAL file');
 
561
        goto 666;
 
562
     end;
 
563
     lTemp4D := CreateDecompressed4D(lImageNamesAll);
 
564
     if (lnSubjAll < 1) or (lnFactors < 1) then begin
 
565
        Showmessage('AnaCOM error: not enough patients ('+inttostr(lnSubjAll)+') or factors ('+inttostr(lnFactors)+').');
 
566
        goto 666;
 
567
     end;
 
568
     lMaskname := lImageNamesAll[0];
 
569
     if not NIFTIhdr_LoadHdr(lMaskname,lMaskHdr) then begin
 
570
           showmessage('Error reading 1st file: '+lMaskName);
 
571
           goto 666;
 
572
     end;
 
573
     lMaskVoxels := ComputeImageDataBytes8bpp(lMaskHdr);
 
574
     if (lMaskVoxels < 2) or (not CheckVoxels(lMaskname,lMaskVoxels,0)){make sure there is uncompressed .img file}  then begin
 
575
           showmessage('Mask file size too small.');
 
576
           goto 666;
 
577
     end;
 
578
     lOutName := ExtractFileDirWithPathDelim(lMaskName)+'results';
 
579
     MainForm.SaveHdrDlg.Filename := loutname;
 
580
     lOutName := lOutName+'.nii.gz';
 
581
     if not MainForm.SaveHdrName ('Base Statistical Map', lOutName) then exit;
 
582
     for lFact := 1 to lnFactors do begin
 
583
         lImageNames.clear;
 
584
         for lSubj := 1 to lnSubjAll do
 
585
           {$IFNDEF FPC}if lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)] <> NaN then {$ENDIF}
 
586
              lImageNames.Add(lImageNamesAll[lSubj-1]);
 
587
         lnSubj := lImageNames.Count;
 
588
         if lnSubj > 1 then begin
 
589
            getmem(lSymptomRA,lnSubj * sizeof(single));
 
590
            lnSubj := 0;
 
591
            for lSubj := 1 to lnSubjAll do
 
592
              {$IFNDEF FPC}if lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)] <> NaN then begin
 
593
              {$ELSE} begin{$ENDIF}
 
594
                 inc(lnSubj);
 
595
                 lSymptomRA^[lnSubj] := lMultiSymptomRA^[lSubj+((lFact-1)*lnSubjAll)];
 
596
              end;
 
597
            MainForm.NPMmsgClear;
 
598
            MainForm.NPMMsg(MainForm.GetKVers);
 
599
            MainForm.NPMMsg('Threads: '+inttostr(gnCPUThreads));
 
600
            npmform.MainForm.memo1.lines.add('Control (n='+inttostr(lnControlObservations)+')performance   ['+lControlFilename+']');
 
601
            for lI := 1 to lnControlObservations do
 
602
                npmform.MainForm.memo1.lines.add(inttostr(lI)+' '+floattostr(lControldata^[lI]));
 
603
            lFactName := lPredictorList.Strings[lFact-1];
 
604
            lFactName := LegitFilename(lFactName,lFact);
 
605
            MainForm.NPMMsg('Patient performance, (n= '+inttostr(lnSubj)+') Factor = '+lFactname);
 
606
            For lSubj := 1 to lnSubj do
 
607
                MainForm.NPMMsg (lImageNames.Strings[lSubj-1] + ' = '+realtostr(lSymptomRA^[lSubj],2) );
 
608
            MainForm.NPMMsg('Total voxels = '+inttostr(lMaskVoxels));
 
609
            MainForm.NPMMsg('Only testing voxels damaged in at least '+inttostr(lnCrit)+' individual[s]');
 
610
            MainForm.NPMMsg('Number of Lesion maps = '+inttostr(lnSubj));
 
611
            if not CheckVoxelsGroup(lImageNames,lMaskVoxels) then begin
 
612
               showmessage('File dimensions differ from mask.');
 
613
               goto 666;
 
614
            end;
 
615
            MainForm.ReportDescriptives(lSymptomRA,lnSubj);
 
616
            AnacomLesionNPMAnalyze(lImageNames,lMaskHdr,lnCrit,-1,lnControlObservations,lSymptomRA,lControldata,lFactName,lOutname,true {ttest},false{BM});
 
617
            Freemem(lSymptomRA);
 
618
         end; //lnsubj > 1
 
619
     end; //for each factor
 
620
     if lnSubjAll > 0 then
 
621
       Freemem(lMultiSymptomRA);
 
622
     666:
 
623
     lImageNames.Free;
 
624
     lImageNamesAll.Free;
 
625
     lPredictorList.Free;
 
626
     DeleteDecompressed4D(lTemp4D);
 
627
     freemem(lControldata);
 
628
end;
 
629
 
 
630
end.