1 /*****************************************************************************
2
3
4 ****************************************************************************/
5
6 #include <stdio.h>
7 #include <dirent.h>
8 #include <malloc.h>
9 #include <string.h>
10 #include <unistd.h>
11 #include <stdbool.h>
12 #include "mrilib.h"
13 #include "distanceField.h"
14 #include <float.h>
15
16 #define PROCESS_ROIS_SEPARATELY 0 // Applies to marching parabolas
17 #define BIG 10000000000.0
18
19 typedef float flt;
20
21 typedef enum METRIC_TYPE {
22 MARCHING_PARABOLAS,
23 EROSION
24 } METRIC_TYPE ;
25
26
27 int Cleanup(char *inputFileName, char *outputFileName, THD_3dim_dataset *din);
28 int afni_edt(THD_3dim_dataset * din, float *outImg,
29 bool do_sqrt, bool edges_are_zero_for_nz, bool debugMode);
30 int erosion(THD_3dim_dataset * din, float *outImg);
31 int open_input_dset(THD_3dim_dataset ** din, char * fname);
32 int outputDistanceField(THD_3dim_dataset *dout, char *outputFileName);
33 int outputDistanceFieldDebug(float *outImg, THD_3dim_dataset *din,
34 char *outputFileName);
35 int doesFileExist(char * searchPath, char * prefix,
36 char *appendage , char * outputFileName);
37 void edt1_local(THD_3dim_dataset * din, flt * df, int n);
38 void edt_local(float scale, flt * f, int n);
39 float sqr(float x);
40 static flt vx(flt * f, int p, int q);
41 bool sixConnectedAllHi(BYTE *buffer, int index, int nx, int ny, int nz);
42 int getIndex(int x, int y, int z, int nx, int ny, int nz);
43 int transposeYZ(float *volume, int nx, int ny, int nz);
44 int testTransposeFunction(THD_3dim_dataset * din);
45 int outputTransposedDataSet(float *buffer, THD_3dim_dataset *din,
46 int nx, int nz, int ny);
47 int shortToByte(THD_3dim_dataset ** din);
48 ERROR_NUMBER getNonzeroIndices(int nvox, int *inputImg,
49 int *numIndices, int **indices);
50 ERROR_NUMBER processIndex(int index, int *inputImg,
51 float **outImg, THD_3dim_dataset *din);
52 int usage();
53 ERROR_NUMBER img3d_Euclidean_DT(int *im, int nx, int ny, int nz,
54 bool do_sqrt, bool edges_are_zero_for_nz,
55 float *ad3, float *odt);
56 ERROR_NUMBER run_EDTD_per_line(int *roi_line, float *dist2_line, int Na,
57 float delta, bool edges_are_zero_for_nz);
58 float * Euclidean_DT_delta(float *f, int n, float delta);
59 ERROR_NUMBER getDistanceFieldDataSet(THD_3dim_dataset *din,
60 THD_3dim_dataset **dout, int metric,
61 bool do_sqrt, bool edges_are_zero_for_nz,
62 bool debugMode);
63
64 // Debugging variables
65 int debugNx, debugNy, debugNz;
66 float debugScaleFactors[3], *debugOutImage;
67
sqr(float x)68 float sqr(float x){
69 return x*x;
70 }
71
main(int argc,char * argv[])72 int main( int argc, char *argv[] )
73 {
74 char *inputFileName=NULL, *outputFileName=NULL;
75 int i, metric=MARCHING_PARABOLAS;
76 THD_3dim_dataset *din = NULL, *dout = NULL;
77 ERROR_NUMBER errorNumber;
78 float *outImg;
79 bool do_sqrt=TRUE, edges_are_zero_for_nz=TRUE, debugMode = FALSE;
80
81 for (i=0; i<argc; ++i) if (argv[i][0]=='-'){
82 switch(argv[i][1]){
83 case 'd': debugMode = TRUE;
84 break;
85 case 'e':
86 edges_are_zero_for_nz = atoi(argv[++i]);
87 break;
88 case 'i':
89 if (!(inputFileName=(char*)malloc(strlen(argv[++i])+8)))
90 return Cleanup(inputFileName, outputFileName, din);
91 sprintf(inputFileName,"%s",argv[i]);
92 break;
93
94 case 'm':
95 if (!strcmp(argv[++i],"MARCHING_PARABOLAS")){
96 metric = MARCHING_PARABOLAS;
97 }
98 if (!strcmp(argv[i],"EROSION")){
99 metric = EROSION;
100 } else {
101 fprintf(stderr, "Unrecognized metric");
102 if (inputFileName) free(inputFileName);
103 return ERROR_UNRECOGNIZABLE_SPECIES;
104 }
105 break;
106
107 case 'o':
108 if (!(outputFileName=(char*)malloc(strlen(argv[++i])+8)))
109 return Cleanup(inputFileName, outputFileName, din);
110 sprintf(outputFileName,"%s",argv[i]);
111 break;
112
113 case 's':
114 do_sqrt = atoi(argv[++i]);
115 break;
116 default:
117 return usage();
118 }
119 }
120
121 if (!inputFileName) return usage();
122
123 // Open dataset
124 if ( (errorNumber=open_input_dset(&din, inputFileName))!=ERROR_NONE ){
125 Cleanup(inputFileName, outputFileName, din);
126 return errorNumber;
127 }
128
129 if ( (errorNumber=getDistanceFieldDataSet(din, &dout, metric, do_sqrt,
130 edges_are_zero_for_nz,
131 debugMode)) != ERROR_NONE ){
132 Cleanup(inputFileName, outputFileName, din);
133 return errorNumber;
134 }
135
136 if (!outputFileName){ // Default output filename
137 char *prefix=DSET_PREFIX(din);
138 char *searchPath=DSET_DIRNAME(din);
139 char appendage[256];
140
141 // Set appendage
142 switch (metric){
143 case MARCHING_PARABOLAS:
144 sprintf(appendage, "MarParab");
145 break;
146 case EROSION:
147 sprintf(appendage, "Erosion");
148 break;
149 }
150
151 if (debugMode) sprintf(appendage, "%s", strcat(appendage, "Debug"));
152
153 // Allocate memory to output name buffer
154 if (!(outputFileName=(char *)malloc(strlen(searchPath) +
155 strlen(prefix) +
156 strlen(appendage)+8))){
157 return ERROR_MEMORY_ALLOCATION;
158 }
159
160 sprintf(outputFileName,"%s%s%s",searchPath,prefix,appendage);
161 }
162
163 if ( (errorNumber=outputDistanceField(dout, outputFileName))!=ERROR_NONE ){
164 Cleanup(inputFileName, outputFileName, din);
165 return errorNumber;
166 }
167
168 Cleanup(inputFileName, outputFileName, din);
169
170 return 0;
171 }
172
usage()173 int usage(){
174 fprintf(stderr,
175 "\n"
176 "This program determines depth of voxels in 3D binary objects, using the\n"
177 "(highly) computationally efficient Euclidean distance transform (EDT).\n"
178 "\n"
179 "SYNOPSIS ~1~\n"
180 "\n"
181 " distanceField -i <input filename> [-o <output filename>][-m <metric>]\n"
182 "\n"
183 "The input file is expected to be a volumetric dataset with integer-valued\n"
184 "voxels.\n"
185 "\n"
186 "The 'metric' is a text string (upper case) describing the algorithm used to\n"
187 "estimate the depth. It may be one of the following.\n"
188 "MARCHING_PARABOLAS - Marching parabolas (default)\n"
189 "EROSION - Erosion algorithm.\n"
190 "\n"
191 "Optional arguments specifically for MARCHING_PARABOLAS:\n"
192 " s: Square root the output\n"
193 " e: Treat edge of field of view as zero (default)\n"
194 " d: (Debug mode.) Generate test object set internally\n"
195 "\n"
196 "If the output filename is not specified, a default name is assigned.\n"
197 "The default name reflects the input filename, metric and whether \n"
198 "the program was run in debug mode.\n"
199 );
200
201 return 0;
202 }
203
getDistanceFieldDataSet(THD_3dim_dataset * din,THD_3dim_dataset ** dout,int metric,bool do_sqrt,bool edges_are_zero_for_nz,bool debugMode)204 ERROR_NUMBER getDistanceFieldDataSet(THD_3dim_dataset *din,
205 THD_3dim_dataset **dout, int metric,
206 bool do_sqrt, bool edges_are_zero_for_nz,
207 bool debugMode){
208
209 ERROR_NUMBER errorNumber;
210 float *outImg;
211
212 if (!(outImg = (float *)calloc(DSET_NVOX(din), sizeof(float)))){
213 return ERROR_MEMORY_ALLOCATION;
214 }
215
216 // Apply metric
217 switch (metric){
218 case MARCHING_PARABOLAS:
219 if ((errorNumber=afni_edt(din, outImg, do_sqrt, edges_are_zero_for_nz,
220 debugMode))!=ERROR_NONE){
221 return errorNumber;
222 }
223 break;
224 case EROSION:
225 if ((errorNumber=erosion(din, outImg))!=ERROR_NONE){
226 return errorNumber;
227 }
228 break;
229 }
230
231 *dout = EDIT_empty_copy(din);
232 EDIT_dset_items( *dout ,
233 ADN_type, MRI_float,
234 ADN_none ) ;
235
236 if (debugMode){
237 THD_ivec3 nxyz={20, 40, 80};
238 THD_fvec3 xyzdel = {2.0, 1.0, 0.5};
239
240 EDIT_dset_items( *dout ,
241 ADN_nxyz, nxyz,
242 ADN_xyzdel, xyzdel,
243 ADN_none ) ;
244 }
245
246 EDIT_substitute_brick(*dout, 0, MRI_float, outImg);
247
248 return ERROR_NONE;
249 }
250
outputDistanceField(THD_3dim_dataset * dout,char * outputFileName)251 int outputDistanceField(THD_3dim_dataset *dout, char *outputFileName){
252
253 // Output Fourier distance image
254 EDIT_dset_items( dout ,
255 ADN_prefix, outputFileName,
256 ADN_none ) ;
257 DSET_write(dout);
258
259 return ERROR_NONE;
260 }
261
outputDistanceFieldDebug(float * outImg,THD_3dim_dataset * din,char * outputFileName)262 int outputDistanceFieldDebug(float *outImg, THD_3dim_dataset *din,
263 char *outputFileName){
264
265 // Set output dimensions
266 THD_ivec3 nxyz={debugNx, debugNy, debugNz};
267 THD_fvec3 xyzdel = {debugScaleFactors[2], debugScaleFactors[1],
268 debugScaleFactors[0]};
269
270 // Output Fourier distance image
271 THD_3dim_dataset *dout = EDIT_empty_copy(din);
272 EDIT_dset_items( dout ,
273 ADN_prefix, outputFileName,
274 ADN_type, MRI_float,
275 ADN_nxyz, nxyz,
276 ADN_xyzdel, xyzdel,
277 ADN_none ) ;
278 EDIT_substitute_brick(dout, 0, MRI_float, outImg);
279 DSET_write(dout);
280
281 return ERROR_NONE;
282 }
283
erosion(THD_3dim_dataset * din,float * outImg)284 int erosion(THD_3dim_dataset * din, float *outImg){
285
286 // Get dimensions in voxels
287 int nz = DSET_NZ(din);
288 int ny = DSET_NY(din);
289 int nx = DSET_NX(din);
290 int nvox = nx*ny*nz;
291 int i;
292 bool objectVoxelsLeft=TRUE;
293 BYTE * buffer;
294
295 if ((nvox < 1) || (nx < 2) || (ny < 2) || (nz < 1))
296 return ERROR_DIFFERENT_DIMENSIONS;
297
298 int brickType=DSET_BRICK_TYPE(din, 0);
299 if( brickType != MRI_byte ) return ERROR_DATATYPENOTHANDLED;
300 BYTE * img = DSET_ARRAY(din, 0);
301
302 // Add uneroded volume to output
303 for (i=0; i<nvox; ++i) outImg[i]+=(img[i]!=0);
304
305 // Allocate memory to buffer
306 if (!(buffer = (BYTE *)malloc(nvox*sizeof(BYTE))))
307 return ERROR_MEMORY_ALLOCATION;
308
309 // Erode volume, adding erode volume to output until no object
310 // voxels left
311 do {
312 objectVoxelsLeft = FALSE;
313
314 // Copy inpout image to buffer
315 for (i=0; i<nvox; ++i) buffer[i] = img[i];
316
317 // Erode input
318 for (i=0; i<nvox; ++i){
319 img[i]=(buffer[i]!=0) && sixConnectedAllHi(buffer, i, nx, ny, nz);
320 if (img[i]) objectVoxelsLeft = TRUE;
321 }
322
323 // Add eroded volume to output
324 for (i=0; i<nvox; ++i) outImg[i]+=(img[i]!=0);
325 } while (objectVoxelsLeft);
326
327 // Cleanup
328 free(buffer);
329
330
331 return ERROR_NONE;
332 }
333
sixConnectedAllHi(BYTE * buffer,int index,int nx,int ny,int nz)334 bool sixConnectedAllHi(BYTE *buffer, int index, int nx, int ny, int nz){
335 int planeSize = nx*ny;
336 int z = (int)(index/planeSize);
337 int y = (int)((index-z*planeSize)/nx);
338 int x = index - z*planeSize - y*nx;
339
340 return buffer[getIndex((x>0)? x-1:x,y,z, nx, ny,nz)] && \
341 buffer[getIndex((x<nx-1)? x+1:x,y,z, nx, ny,nz)] && \
342 buffer[getIndex(x, (y>0)? y-1:y,z, nx, ny,nz)] && \
343 buffer[getIndex(x, (y<ny-1)? y+1:y,z, nx, ny,nz)] && \
344 buffer[getIndex(x,y,(z>0)? z-1:z, nx, ny,nz)] && \
345 buffer[getIndex(x,y,(z<nz-1)? z+1:z, nx, ny,nz)];
346 }
347
getIndex(int x,int y,int z,int nx,int ny,int nz)348 int getIndex(int x, int y, int z, int nx, int ny, int nz){
349 return z*nx*ny + y*nx + x;
350 }
351
352
afni_edt(THD_3dim_dataset * din,float * outImg,bool do_sqrt,bool edges_are_zero_for_nz,bool debugMode)353 int afni_edt(THD_3dim_dataset * din, float *outImg, bool do_sqrt,
354 bool edges_are_zero_for_nz,
355 bool debugMode){
356
357 // Get dimensions in voxels
358 int nz = DSET_NZ(din);
359 int ny = DSET_NY(din);
360 int nx = DSET_NX(din);
361 int nvox = nx*ny*nz;
362 int *inputImg;
363 int *indices;
364 BYTE * byteImg;
365 short * shortImg;
366 float *floatImg, ad3[3];
367 int *vol;
368
369 if ((nvox < 1) || (nx < 2) || (ny < 2) || (nz < 1))
370 return ERROR_DIFFERENT_DIMENSIONS;
371
372 // Alliocate memory to integer input buffer
373 if (!(inputImg=(int *)calloc(nvox,sizeof(int))))
374 return ERROR_MEMORY_ALLOCATION;
375
376 DSET_load(din);
377 int brickType=DSET_BRICK_TYPE(din, 0);
378 fprintf(stderr, "brickType=%d\n", brickType);
379 switch(brickType){
380 case MRI_byte:
381 byteImg = DSET_ARRAY(din, 0);
382 for (int i=0; i<nvox; ++i) inputImg[i] = byteImg[i];
383 break;
384 case MRI_short:
385 shortImg = DSET_ARRAY(din, 0);
386 for (int i=0; i<nvox; ++i) inputImg[i] = shortImg[i];
387 break;
388 case MRI_int:
389 free(inputImg);
390 inputImg = DSET_ARRAY(din, 0);
391 break;
392 case MRI_float:
393 floatImg = DSET_ARRAY(din, 0);
394 for (int i=0; i<nvox; ++i) inputImg[i] = (int)(floatImg[i]+0.5);
395 break;
396 }
397
398 if (debugMode){
399 // DEBUG: Make test volume
400 // make a test image: a map of various ROIs
401 debugNz = nz = 80;
402 debugNy = ny = 40;
403 debugNx = nx = 20;
404 nvox = nx*ny*nz;
405 vol = (int *)calloc(nx*ny*nz, sizeof(int));
406 //free(outImg);
407 outImg = (float *)calloc(nvox,sizeof(float));
408 debugOutImage = outImg;
409 ad3[0] = 0.5;
410 ad3[1] = 1.0;
411 ad3[2] = 2.0;
412 memcpy(debugScaleFactors, ad3, 3*sizeof(float));
413
414 // LX, LY, LZ = np.shape(vol)
415
416 int planesize=nx*ny;
417
418 for (int z=4; z<18; ++z){
419 int zOffset=z*planesize;
420 for (int y=2; y<8; ++y){
421 int yOffset = zOffset + y*nx;
422 for (int x=4; x<8; ++x){
423 vol[yOffset+x] = 1;
424 }
425 }
426 }
427
428 for (int z=2; z<8; ++z){
429 int zOffset=z*planesize;
430 for (int y=4; y<18; ++y){
431 int yOffset = zOffset + y*nx;
432 for (int x=4; x<8; ++x){
433 vol[yOffset+x] = 4;
434 }
435 }
436 }
437
438 for (int z=21; z<30; ++z){
439 int zOffset=z*planesize;
440 for (int y=0; y<10; ++y){
441 int yOffset = zOffset + y*nx;
442 for (int x=11; x<14; ++x){
443 vol[yOffset+x] = 7;
444 }
445 }
446 }
447
448 for (int z=0; z<nz; ++z){
449 int zOffset=z*planesize;
450 for (int y=0; y<ny; ++y){
451 int yOffset = zOffset + y*nx;
452 for (int x=0; x<nx; ++x){
453 if (pow((15-x),2.0) + pow((25-y),2.0) + pow((10-z),2.0)< 31)
454 vol[yOffset+x] = 1;
455 }
456 }
457 }
458
459 for (int z=17; z<19; ++z){
460 int zOffset=z*planesize;
461 for (int y=0; y<ny; ++y){
462 int yOffset = zOffset + y*nx;
463 for (int x=3; x<6; ++x){
464 vol[yOffset+x] = 1;
465 }
466 }
467 }
468
469 for (int z=0; z<nz; ++z){
470 int zOffset=z*planesize;
471 for (int y=19; y<21; ++y){
472 int yOffset = zOffset + y*nx;
473 for (int x=3; x<6; ++x){
474 vol[yOffset+x] = 10;
475 }
476 }
477 }
478
479 for (int z=60; z<70; ++z){
480 int zOffset=z*planesize;
481 for (int y=28; y<36; ++y){
482 int yOffset = zOffset + y*nx;
483 for (int x=14; x<19; ++x){
484 vol[yOffset+x] = 17;
485 }
486 }
487 }
488
489 for (int z=0; z<nz; ++z){
490 int zOffset=z*planesize;
491 for (int y=0; y<ny; ++y){
492 int yOffset = zOffset + y*nx;
493 for (int x=0; x<nx; ++x){
494 if ((pow((65-x),2.0) + pow((10-y),2.0) + pow((10-z),2.0)< 75) &&
495 !(pow((60-x),2.0) + pow((7-y),2.0) + pow((10-z),2.0)< 31))
496 vol[yOffset+x] = 2;
497 }
498 }
499 }
500
501 for (int z=40; z<56; ++z){
502 int zOffset=z*planesize;
503 for (int y=25; y<33; ++y){
504 int yOffset = zOffset + y*nx;
505 for (int x=4; x<8; ++x){
506 vol[yOffset+x] = 5;
507 }
508 }
509 }
510 }
511
512 #if PROCESS_ROIS_SEPARATELY
513 // Get unique nonzero index values
514 if ( (errorNumber=getNonzeroIndices(nvox, inputImg, &numIndices, &indices)) \
515 != ERROR_NONE ){
516 free(inputImg);
517 return errorNumber;
518 }
519
520 // Process each index
521 for (int i=0; i<numIndices; ++i){
522 if ( (errorNumber=processIndex(indices[i], inputImg, &addend, din)) \
523 != ERROR_NONE ){
524 free(inputImg);
525 free(indices);
526 return errorNumber;
527 }
528 for (int j=0; j<nvox; ++j) outImg[j]+=addend[j];
529 }
530
531 free(indices);
532 #else
533 if (debugMode){
534 inputImg = vol;
535 } else {
536 // Get real world voxel sizes
537 ad3[0] = fabs(DSET_DX(din));
538 ad3[1] = fabs(DSET_DY(din));
539 ad3[2] = fabs(DSET_DZ(din));
540 }
541
542 img3d_Euclidean_DT(inputImg, nx, ny, nz,
543 do_sqrt, edges_are_zero_for_nz, ad3, outImg);
544 #endif
545
546 // Cleanup
547 free(inputImg);
548
549 return ERROR_NONE;
550 }
551
img3d_Euclidean_DT(int * im,int nx,int ny,int nz,bool do_sqrt,bool edges_are_zero_for_nz,float * ad3,float * odt)552 ERROR_NUMBER img3d_Euclidean_DT(int *im, int nx, int ny, int nz,
553 bool do_sqrt, bool edges_are_zero_for_nz,
554 float *ad3, float *odt){
555
556 int i, x, y, z;
557 int offset, planeOffset;
558 int nvox = nx*ny*nz; // number of voxels
559 int planeSize = nx*ny;
560
561 // initialize the "output" or answer array
562 // for (int i=0; i<nvox; ++i) odt[i] = (im[i]>0)? BIG : 0;
563 for ( i=0; i<nvox; ++i )
564 odt[i] = BIG;
565
566 // first pass: start with all BIGs (along x-axes)
567 int *inRow = im;
568 float *outRow = odt;
569 for ( z = 0; z <nz; ++z ){
570 for ( y = 0; y < ny; ++y ){
571 // Calc with it, and save results
572 run_EDTD_per_line( inRow, outRow, nx, ad3[0],
573 edges_are_zero_for_nz );
574
575 // Increment row
576 inRow += nx;
577 outRow += nx;
578 }
579 }
580
581 if (!(inRow=(int *)malloc(ny*sizeof(int))))
582 return ERROR_MEMORY_ALLOCATION;
583 if (!(outRow=(float *)calloc(ny,sizeof(float)))) {
584 free(inRow);
585 return ERROR_MEMORY_ALLOCATION;
586 }
587 for ( z = 0; z < nz; ++z ){
588 int *inPlane = im + (z*planeSize);
589 float *outPlane = odt + (z*planeSize);
590 for ( x = 0; x < nx; ++x ){
591 // get a line...
592 for ( y=0; y<ny; ++y ){
593 offset = (y*nx) + x;
594 inRow[y] = inPlane[offset ];
595 outRow[y] = outPlane[offset ];
596 }
597
598 // ... and then calc with it, and save results
599 run_EDTD_per_line( inRow, outRow, ny, ad3[1],
600 edges_are_zero_for_nz );
601
602 // Record new output row
603 for ( y=0; y<ny; ++y ){
604 offset = (y*nx) + x;
605 outPlane[offset] = outRow[y];
606 }
607 }
608 }
609 free(outRow);
610 free(inRow);
611
612 // 2nd pass: start from previous; any other dimensions would carry
613 // on from here
614 if (!(inRow=(int *)malloc(nz*sizeof(int))))
615 return ERROR_MEMORY_ALLOCATION;
616 if (!(outRow=(float *)calloc(nz, sizeof(float)))) {
617 free(inRow);
618 return ERROR_MEMORY_ALLOCATION;
619 }
620 for ( y = 0; y < ny; ++y ) {
621 for ( x = 0; x < nx; ++x ){
622 planeOffset = (y*nx) + x;
623 // get a line...
624 for (int z=0; z<nz; ++z){
625 offset = planeOffset + (z*planeSize);
626 inRow[z] = im[offset] ;
627 outRow[z] = odt[offset] ;
628 }
629
630 // ... and then calc with it, and save results
631 run_EDTD_per_line( inRow, outRow, nz, ad3[2],
632 edges_are_zero_for_nz );
633
634 // Record new output row
635 for ( z=0; z<nz; ++z ) {
636 offset = planeOffset + (z*planeSize);
637 odt[offset] = outRow[z];
638 }
639 }
640 }
641 free(outRow);
642 free(inRow);
643
644 if (do_sqrt)
645 for ( i=0; i<nvox; ++i )
646 odt[i] = sqrt(odt[i]);
647
648 return ERROR_NONE;
649 }
650
run_EDTD_per_line(int * roi_line,float * dist2_line,int Na,float delta,bool edges_are_zero_for_nz)651 ERROR_NUMBER run_EDTD_per_line(int *roi_line, float *dist2_line, int Na,
652 float delta, bool edges_are_zero_for_nz) {
653 int idx = 0;
654 int i, m, n;
655 float *line_out=NULL;
656 int start, stop, inc, roi;
657
658 float *Df = NULL;
659
660 int limit = Na-1;
661 size_t rowLengthInBytes = Na*sizeof(float);
662
663 if (!(line_out=(float *)malloc(rowLengthInBytes)))
664 return ERROR_MEMORY_ALLOCATION;
665
666 while (idx < Na){
667 // get interval of line with current ROI value
668 roi = roi_line[idx];
669 n = idx;
670 while (n < Na){
671 if (roi_line[n] != roi){
672 break;
673 }
674 n += 1;
675 }
676 n -= 1;
677 // n now has the index of last matching element
678
679 float *paddedLine=(float *)calloc(Na+2,sizeof(float));
680 // actual ROI is in range of indices [start, stop] in
681 // paddedLine. 'inc' will tell us length of distance array
682 // put into Euclidean_DT_delta(), which can include padding at
683 // either end.
684 start = 0;
685 stop = limit;
686 inc = 0;
687
688 if (idx != 0 || (edges_are_zero_for_nz && roi != 0)){
689 start = 1;
690 inc = 1;
691 }
692 // put actual values from dist**2 field...
693 for ( m=idx; m<=n; ++m ){
694 paddedLine[inc++] = dist2_line[m];
695 }
696 // inc finishes 1 greater than the actual end: is length so far
697 stop = inc-1; // 'stop' is index of actual end of roi
698 // pad at end?
699 if (n < limit || (edges_are_zero_for_nz && roi != 0)){
700 inc+=1;
701 }
702
703 // [PT] and 'inc' should already have correct value from
704 // above; don't add 1 here
705 Df = Euclidean_DT_delta(paddedLine, inc, delta);
706
707 memcpy(&(line_out[idx]), &(Df[start]), (stop-start+1)*sizeof(float));
708
709 free(paddedLine);
710 if( Df )
711 free(Df);
712
713 idx = n+1;
714 }
715
716 // DEBUG
717 for ( i=0; i<Na; ++i )
718 if (line_out[i]==0){
719 fprintf(stderr, "Zero valued distance\n");
720 }
721
722 memcpy(dist2_line, line_out, rowLengthInBytes);
723 free(line_out);
724
725 return ERROR_NONE;
726 }
727
Euclidean_DT_delta(float * f,int n,float delta)728 float * Euclidean_DT_delta(float *f, int n, float delta){
729 // Classical Euclidean Distance Transform (EDT) of Felzenszwalb and
730 // Huttenlocher (2012), but for given voxel lengths.
731 //
732 // Assumes that: len(f) < sqrt(10**10).
733 //
734 // In this version, all voxels should have equal length, and units
735 // are "edge length" or "number of voxels."
736 //
737 // Parameters
738 // ----------
739 //
740 // f : 1D array or list, distance**2 values (or, to start, binarized
741 // between 0 and BIG).
742 //
743 // delta : voxel edge length size along a particular direction
744
745 int q;
746 int *v=NULL;
747 int k = 0;
748 float *z=NULL, *Df=NULL;
749 float s;
750
751 if (!(v=(int *)calloc(n, sizeof(int))) ||
752 !(z=(float *)calloc(n+1, sizeof(float)))){
753 if (v) free(v);
754 return NULL;
755 }
756 z[0] = -BIG;
757 z[1] = BIG;
758
759 for ( q = 1; q<n; ++q ) {
760 s = f[q] + pow(q*delta, 2.0) - (f[v[k]] + pow(v[k]*delta,2.0));
761 s/= 2. * delta * (q - v[k]);
762 while (s <= z[k]){
763 k-= 1;
764 s = f[q] + pow(q*delta,2.0) - (f[v[k]] + pow(v[k]*delta, 2.0));
765 s/= 2. * delta * (q - v[k]);
766 }
767 k+= 1;
768 v[k] = q;
769 z[k] = s;
770 z[k+1] = BIG;
771 }
772
773 k = 0;
774 if (!(Df=(float *)calloc(n, sizeof(float)))){
775 free(v);
776 free(z);
777 return NULL;
778 }
779 for ( q=0; q<n; ++q ){
780 while (z[k+1] < q * delta) k+= 1;
781 Df[q] = pow(delta*(q - v[k]), 2.0) + f[v[k]];
782 }
783
784 free(v);
785 free(z);
786
787 return Df;
788 }
789
processIndex(int index,int * inputImg,float ** outImg,THD_3dim_dataset * din)790 ERROR_NUMBER processIndex(int index, int *inputImg, float **outImg,
791 THD_3dim_dataset *din){
792 // Get dimensions in voxels
793 int r, v, x, y, z;
794 int xo, yo, zo;
795 size_t i, vo, nRow;
796 int nz = DSET_NZ(din);
797 int ny = DSET_NY(din);
798 int nx = DSET_NX(din);
799 int nvox = nx*ny*nz;
800
801 int nVol = 1;
802 int nvox3D = nx * ny * MAX(nz, 1);
803
804 nVol = nvox / nvox3D;
805 if ((nvox3D * nVol) != nvox)
806 return ERROR_DIFFERENT_DIMENSIONS;
807
808 // Get real world voxel sizes
809 float xDim = fabs(DSET_DX(din));
810 float yDim = fabs(DSET_DY(din));
811 float zDim = fabs(DSET_DZ(din));
812
813 float yDimSqrd = yDim*yDim;
814 float zDimSqrd = zDim*zDim;
815
816 if (!(*outImg=(float *)calloc(nvox, sizeof(float))))
817 return ERROR_MEMORY_ALLOCATION;
818
819 for ( i = 0; i < nvox; i++ ) {
820 if (inputImg[i] == index)
821 (*outImg)[i] = BIG;
822 }
823 nRow = ny*nz;
824
825 //EDT in left-right direction
826 for ( r = 0; r < nRow; r++ ) {
827 flt * imgRow = (*outImg) + (r * nx);
828 // edt1_local(din, imgRow, nx);
829 edt_local(xDim, imgRow, nx);
830 }
831
832 //EDT in anterior-posterior direction
833 nRow = nx * nz; //transpose XYZ to YXZ and blur Y columns with XZ Rows
834 for ( v = 0; v < nVol; v++ ) { //transpose each volume separately
835 //alloc for each volume to allow openmp
836 flt * img3D = (flt *)calloc(nvox3D*sizeof(flt), 64);
837
838 //transpose data
839 vo = v * nvox3D; //volume offset
840 for ( z = 0; z < nz; z++ ) {
841 zo = z * nx * ny;
842 for ( y = 0; y < ny; y++ ) {
843 xo = 0;
844 for ( x = 0; x < nx; x++ ) {
845 img3D[zo+xo+y] = (*outImg)[vo]/yDimSqrd;
846 vo += 1;
847 xo += ny;
848 }
849 }
850 }
851 //perform EDT for all rows
852 for ( r = 0; r < nRow; r++ ) {
853 flt * imgRow = img3D + (r * ny);
854 edt_local(yDim, imgRow, ny);
855 }
856 //transpose data back
857 vo = v * nvox3D; //volume offset
858 for ( z = 0; z < nz; z++ ) {
859 zo = z * nx * ny;
860 for ( y = 0; y < ny; y++ ) {
861 xo = 0;
862 for ( x = 0; x < nx; x++ ) {
863 (*outImg)[vo] = img3D[zo+xo+y];
864 vo += 1;
865 xo += ny;
866 }
867 }
868 }
869 free (img3D);
870 } //for each volume
871
872 //EDT in head-foot direction
873 nRow = nx * ny; //transpose XYZ to ZXY and blur Z columns with XY Rows
874 for ( v = 0; v < nVol; v++ ) { //transpose each volume separately
875 //alloc for each volume to allow openmp
876 flt * img3D = (flt *)calloc(nvox3D*sizeof(flt), 64);
877 //transpose data
878 vo = v * nvox3D; //volume offset
879 for ( z = 0; z < nz; z++ ) {
880 for ( y = 0; y < ny; y++ ) {
881 yo = y * nz * nx;
882 xo = 0;
883 for ( x = 0; x < nx; x++ ) {
884 img3D[z+xo+yo] = (*outImg)[vo]/zDimSqrd;
885 vo += 1;
886 xo += nz;
887 }
888 }
889 }
890 //perform EDT for all "rows"
891 for ( r = 0; r < nRow; r++ ) {
892 flt * imgRow = img3D + (r * nz);
893 edt_local(zDim, imgRow, nz);
894 }
895 //transpose data back
896 vo = v * nvox3D; //volume offset
897 for ( z = 0; z < nz; z++ ) {
898 for ( y = 0; y < ny; y++ ) {
899 yo = y * nz * nx;
900 xo = 0;
901 for ( x = 0; x < nx; x++ ) {
902 (*outImg)[vo] = img3D[z+xo+yo];
903 vo += 1;
904 xo += nz;
905 } //x
906 } //y
907 } //z
908 free (img3D);
909 } //for each volume
910
911 return ERROR_NONE;
912 }
913
getNonzeroIndices(int nvox,int * inputImg,int * numIndices,int ** indices)914 ERROR_NUMBER getNonzeroIndices(int nvox, int *inputImg, int *numIndices,
915 int **indices){
916 int *buffer;
917 int i, j, voxelValue;
918 bool old;
919
920 // Initialize
921 *numIndices = 0;
922 if (!(buffer=(int *)calloc(nvox,sizeof(int))))
923 return ERROR_MEMORY_ALLOCATION;
924
925 // Count unique indices and fill buffer with set of indices
926 for (i=0; i<nvox; ++i) if ((voxelValue=inputImg[i])>0){
927 old = false;
928 // for (j=0; j<*numIndices; ++j)
929 // if ((*indices)[j]==voxelValue)
930 // old = true;
931 // Core dump
932 for (j=0; j<*numIndices; ++j) if (buffer[j]==voxelValue) old = true;
933 if (!old){
934 buffer[(*numIndices)++]=voxelValue;
935 }
936 }
937
938 // Trim output buffer to number of indices
939 if (!(*indices=(int *)malloc(*numIndices*sizeof(int)))){
940 free(buffer);
941 return ERROR_MEMORY_ALLOCATION;
942 }
943 for (i=0; i<*numIndices; ++i) (*indices)[i] = buffer[i];
944 free(buffer);
945
946 return ERROR_NONE;
947 }
948
transposeYZ(float * volume,int nx,int ny,int nz)949 int transposeYZ(float *volume, int nx, int ny, int nz){
950 float * buffer;
951 int nvox=nx*ny*nz;
952
953 if (!(buffer = (float *)malloc(nvox*sizeof(float)))!=ERROR_NONE)
954 return ERROR_MEMORY_ALLOCATION;
955
956 int z0 = 0;
957 for (int z=0; z<nz; ++z){
958 int Y = z0;
959 for (int y=0; y<ny; ++y){
960 for (int x=0; x<nx; ++x){
961 buffer[z0+x]=volume[Y];
962 buffer[Y++]=volume[z0+x];
963 }
964 }
965 z0+=nx*ny;
966 }
967 memcpy((void *)volume, (void *)buffer, nvox*sizeof(float));
968 free(buffer);
969
970 return ERROR_NONE;
971 }
972
testTransposeFunction(THD_3dim_dataset * din)973 int testTransposeFunction(THD_3dim_dataset * din){
974 // Get dimensions in voxels
975 int nz = DSET_NZ(din);
976 int ny = DSET_NY(din);
977 int nx = DSET_NX(din);
978 int nvox = nx*ny*nz;
979 int i, errorNumber;
980 float *buffer;
981
982 // Allocate memory to buffer
983 if (!(buffer = (float *)malloc(nvox*sizeof(float))))
984 return ERROR_MEMORY_ALLOCATION;
985
986 // Read image data
987 BYTE * img = DSET_ARRAY(din, 0);
988
989 for (i=0; i<nvox; ++i) buffer[i]=img[i];
990
991 // Apply transpose
992 if ((errorNumber=transposeYZ(buffer, nx, ny, nz))!=ERROR_NONE){
993 free(buffer);
994 return errorNumber;
995 }
996
997 // Output transpoed volume to afni dataset
998 outputTransposedDataSet(buffer, din, nx, nz, ny);
999
1000 free(buffer);
1001 return ERROR_NONE;
1002 }
1003
1004
outputTransposedDataSet(float * buffer,THD_3dim_dataset * din,int nx,int nz,int ny)1005 int outputTransposedDataSet(float *buffer, THD_3dim_dataset *din,
1006 int nx, int nz, int ny){
1007 char *prefix=DSET_PREFIX(din);
1008 char *searchPath=DSET_DIRNAME(din);
1009 char *outputFileName;
1010 char appendage[]="Transpose";
1011
1012 // Allocate memory to output name buffer
1013 if (!(outputFileName=(char *)malloc(strlen(searchPath)+strlen(prefix)+strlen(appendage)+8))){
1014 return ERROR_MEMORY_ALLOCATION;
1015 }
1016
1017 // Determine whether output file already exists
1018 int outputFileExists = doesFileExist(searchPath,
1019 prefix,
1020 appendage,
1021 outputFileName);
1022
1023 // Set output dimensions
1024 THD_ivec3 nxyz[3]={nx, ny, nz};
1025
1026 // Output Fourier spectrum image (if it does not already exist)
1027 if (!outputFileExists){
1028 THD_3dim_dataset *dout = EDIT_empty_copy(din);
1029 sprintf(outputFileName,"%s%s%s",searchPath,prefix,appendage);
1030 EDIT_dset_items( dout ,
1031 ADN_prefix, outputFileName,
1032 ADN_type, MRI_float,
1033 ADN_nxyz, nxyz,
1034 ADN_none ) ;
1035 EDIT_substitute_brick(dout, 0, MRI_float, buffer);
1036 DSET_write(dout);
1037 }
1038
1039 // Cleanup
1040 free(outputFileName);
1041
1042 return ERROR_NONE;
1043 }
1044
edt1_local(THD_3dim_dataset * din,flt * df,int n)1045 void edt1_local(THD_3dim_dataset * din, flt * df, int n) {
1046 //first dimension is simple
1047
1048 // Get real world voxel sizes
1049 float xDim = fabs(DSET_DX(din));
1050 float yDim = fabs(DSET_DY(din));
1051
1052 xDim=yDim/xDim;
1053 yDim=1.0f;
1054
1055 int q, prevX;
1056 flt prevY, v;
1057 prevX = 0;
1058 prevY = BIG;
1059 //forward
1060 for (q = 0; q < n; q++ ) {
1061 if (df[q] == 0) {
1062 prevX = q;
1063 prevY = 0;
1064 } else
1065 df[q] = sqr((q-prevX)/xDim)+(prevY/yDim);
1066 }
1067 //reverse
1068 prevX = n;
1069 prevY = BIG;
1070 for (q = (n-1); q >= 0; q-- ) {
1071 v = sqr((q-prevX)/xDim)+(prevY/yDim);
1072 if (df[q] < v) {
1073 prevX = q; // DO NOT CHANGE
1074 prevY = df[q]; // DO NOT CHANGE
1075 } else
1076 df[q] = v;
1077 }
1078 }
1079
vx(flt * f,int p,int q)1080 flt vx(flt * f, int p, int q) {
1081 if ((f[p] == BIG) || (f[q] == BIG))
1082 return BIG;
1083 else
1084 return ((f[q] + q*q) - (f[p] + p*p)) / (2.0*q - 2.0*p);
1085 }
1086
edt_local(float scale,flt * f,int n)1087 void edt_local(float scale, flt * f, int n) {
1088 float scaleSqrd = scale*scale;
1089
1090 int q, p, k;
1091 flt s, dx;
1092 flt * d = (flt *)calloc((n)*sizeof(flt), 64);
1093 flt * z = (flt *)calloc((n)*sizeof(flt), 64);
1094 int * v = (int *)calloc((n)*sizeof(int), 64);
1095
1096 // # Find the lower envelope of a sequence of parabolas.
1097 // # f...source data (returns the Y of the parabola vertex at X)
1098 // # d...destination data (final distance values are written here)
1099 // # z...temporary used to store X coords of parabola intersections
1100 // # v...temporary used to store X coords of parabola vertices
1101 // # i...resulting X coords of parabola vertices
1102 // # n...number of pixels in "f" to process
1103 // # Always add the first pixel to the enveloping set since it is
1104 // # obviously lower than all parabolas processed so far.
1105
1106 k = 0;
1107 v[0] = 0;
1108 z[0] = -BIG;
1109 z[1] = BIG;
1110
1111 for (q = 1; q < n; q++ ) {
1112 // If the new parabola is lower than the right-most parabola in
1113 // # the envelope, remove it from the envelope. To make this
1114 // # determination, find the X coordinate of the intersection (s)
1115 // # between the parabolas with vertices at (q,f[q]) and (p,f[p]).
1116 p = v[k]; // DO NOT CHANGE
1117 s = vx(f, p,q);
1118 while (s <= z[k]) {
1119 k = k - 1;
1120 p = v[k]; // DO NOT CHANGE
1121 s = vx(f, p,q); // DO NOT CHANGE
1122 }
1123 //# Add the new parabola to the envelope.
1124 k = k + 1;
1125 v[k] = q; // DO NOT CHANGE
1126 z[k] = s;
1127 z[k + 1] = BIG;
1128 }
1129 // # Go back through the parabolas in the envelope and evaluate them
1130 // # in order to populate the distance values at each X coordinate.
1131 k = 0;
1132 for (q = 0; q < n; q++ ) {
1133 while (z[k + 1] < q)
1134 k = k + 1;
1135 dx = (q - v[k])*scale;
1136 // d[q] = dx * dx + f[v[k]];
1137 d[q] = dx * dx + f[v[k]]*scaleSqrd;
1138 }
1139 for (q = 0; q < n; q++ )
1140 f[q] = d[q];
1141 free (d);
1142 free (z);
1143 free (v);
1144 }
1145
doesFileExist(char * searchPath,char * prefix,char * appendage,char * outputFileName)1146 int doesFileExist(char * searchPath, char * prefix,char *appendage ,
1147 char * outputFileName){
1148 int outputFileExists=0;
1149
1150 struct dirent *dir;
1151 DIR *d;
1152 d = opendir(searchPath);
1153 sprintf(outputFileName,"%s%s",prefix,appendage);
1154 if (d) {
1155 while ((dir = readdir(d)) != NULL) {
1156 if (strstr(dir->d_name,outputFileName)){
1157 outputFileExists=1;
1158 break;
1159 }
1160 }
1161 closedir(d);
1162 }
1163
1164 return outputFileExists;
1165 }
1166
open_input_dset(THD_3dim_dataset ** din,char * fname)1167 int open_input_dset(THD_3dim_dataset ** din, char * fname)
1168 {
1169 *din = THD_open_dataset(fname);
1170 if( ! *din ) {
1171 fprintf(stderr,"** failed to read input dataset '%s'\n", fname);
1172 return ERROR_READING_FILE;
1173 }
1174
1175 return ERROR_NONE;
1176 }
1177
shortToByte(THD_3dim_dataset ** din)1178 int shortToByte(THD_3dim_dataset ** din){
1179 BYTE * outImg;
1180 int i;
1181
1182 // Get dimensions in voxels
1183 int nvox = DSET_NVOX(*din);
1184
1185 if (!(outImg = (BYTE *)calloc(nvox, sizeof(BYTE)))){
1186 return ERROR_MEMORY_ALLOCATION;
1187 }
1188
1189 // Get input data
1190 int brickType=DSET_BRICK_TYPE(*din, 0);
1191 if( brickType != MRI_short ) return ERROR_DATATYPENOTHANDLED;
1192 DSET_load(*din);
1193 short * img = (short *)DSET_ARRAY(*din, 0);
1194
1195 // Map short input into BYTE output
1196 for (i=0; i<nvox; ++i) outImg[i]+=(img[i]!=0);
1197
1198 THD_3dim_dataset *dout = EDIT_empty_copy(*din);
1199 char *prefix=DSET_PREFIX(*din);
1200 char *searchPath=DSET_DIRNAME(*din);
1201 EDIT_dset_items( dout ,
1202 ADN_prefix, strcat(searchPath, prefix),
1203 ADN_type, MRI_byte,
1204 ADN_none ) ;
1205 EDIT_substitute_brick(dout, 0, MRI_byte, outImg);
1206
1207 // Cleanup (data set)
1208 DSET_delete(*din);
1209 *din = dout;
1210
1211 return ERROR_NONE;
1212 }
1213
Cleanup(char * inputFileName,char * outputFileName,THD_3dim_dataset * din)1214 int Cleanup(char *inputFileName, char *outputFileName, THD_3dim_dataset *din){
1215
1216 if (inputFileName) free(inputFileName);
1217 if (outputFileName) free(outputFileName);
1218
1219 return 0;
1220 }
1221