1 #include <stdio.h>
2 #include <stdlib.h>
3 #include <string.h>
4 #include <math.h>
5 #include "mrilib.h"
6 #include "cluster.h"
7
8
9
10 // from command.c file
11
display_help(void)12 static void display_help(void)
13 { printf ("Clustering 4segmentation, command-line version.\n");
14 printf ("USAGE: cluster [options]\n");
15 printf ("options:\n");
16 printf (" -v, --version Version information\n");
17 printf (" -f filename File loading\n");
18 printf (" -cg a|m Specifies whether to center each row\n"
19 " in the data\n"
20 " a: Subtract the mean of each row\n"
21 " m: Subtract the median of each row\n"
22 " (default is no centering)\n");
23 printf (" -ng Specifies to normalize each row in the data\n"
24 " (default is no normalization)\n");
25 printf (" -ca a|m Specifies whether to center each column \n"
26 " in the data\n"
27 " a: Subtract the mean of each column\n"
28 " m: Subtract the median of each column\n"
29 " (default is no centering)\n");
30 printf (" -na Specifies to normalize each column in the data\n"
31 " (default is no normalization)\n");
32 printf (" -u jobname Allows you to specify a different name for the \n"
33 " output files.\n"
34 " (default is derived from the input file name)\n");
35 printf (" -g [0..8] Specifies the distance measure for gene clustering \n"
36 );
37 printf (" 0: No gene clustering\n"
38 " 1: Uncentered correlation\n"
39 " 2: Pearson correlation\n"
40 " 3: Uncentered correlation, absolute value\n"
41 " 4: Pearson correlation, absolute value\n"
42 " 5: Spearman's rank correlation\n"
43 " 6: Kendall's tau\n"
44 " 7: Euclidean distance\n"
45 " 8: City-block distance\n"
46 " (default: 1)\n");
47 printf (" -k number Specifies whether to run k-means clustering\n"
48 " instead of hierarchical clustering, and the number\n"
49 " of clusters k to use\n");
50 printf (" -c number Specifies the number of clusters for tree cutting\n"
51 " after hierarchical clustering\n");
52 printf (" -r number For k-means clustering, the number of times the\n"
53 " k-means clustering algorithm is run\n"
54 " (default: 1)\n");
55 printf (" -m [msca] Specifies which hierarchical clustering method to use\n"
56 " m: Pairwise complete-linkage\n"
57 " s: Pairwise single-linkage\n"
58 " c: Pairwise centroid-linkage\n"
59 " a: Pairwise average-linkage\n"
60 " (default: m)\n");
61 return;
62 }
63
64
display_version(void)65 static void display_version(void)
66 { printf ("\n"
67 "Clustering for segmentation, command line version ,\n"
68 "using the C Clustering Library version 1.38.\n"
69 "\n"
70 "Cluster was originally written by Michael Eisen (eisen 'AT' rana.lbl.gov)\n"
71 "Copyright 1998-99 Stanford University.\n");
72 printf ("\n"
73 "The command line version of Clustering version was created by Andrej Vovk\n"
74 "and will be surely corrected by Ziad Saad ;),\n"
75 "\n");
76 return;
77 }
78
setjobname(const char * basename,int strip)79 static char* setjobname(const char* basename, int strip)
80 { char* jobname;
81 int n = strlen(basename);
82 if (strip)
83 { char* extension = strrchr(basename, '.');
84 if (extension) n -= strlen(extension);
85 }
86 jobname = malloc((n+1)*sizeof(char));
87 strncpy (jobname, basename, n);
88 jobname[n] = '\0';
89 return jobname;
90 }
91
readnumber(const char word[])92 static int readnumber(const char word[])
93 { char* error = 0;
94 long value = 0;
95 if (word[0] == '0' && word[1] != '\0') ++word;
96 value = strtol(word,&error,0);
97 if (*error=='\0') return (int)value;
98 else return -1;
99 }
100
getmetric(int i)101 static char getmetric(int i)
102 { switch (i)
103 { case 1: return 'u';
104 case 2: return 'c';
105 case 3: return 'x';
106 case 4: return 'a';
107 case 5: return 's';
108 case 6: return 'k';
109 case 7: return 'e';
110 case 8: return 'b';
111 default: return '\0';
112 }
113 /* Never get here */
114 return '\0';
115 }
116
117
118
119
120 /* ========================================================================= */
121
122
123
example_distance_gene(int nrows,int ncols,double ** data)124 double** example_distance_gene(int nrows, int ncols, double** data)
125 /* Calculate the distance matrix between genes using the Euclidean distance. */
126 { int i, j, ii, nl, nc;
127 double** distMatrix;
128 double* weight = malloc(ncols*sizeof(double));
129 int** mask = NULL;
130
131 mask = (int **)calloc(sizeof(int*), nrows);
132 for (ii=0;ii<nrows;++ii) {
133 mask[ii] = (int *)calloc(sizeof(int),ncols);
134 }
135
136 for (nl=0; nl<nrows; ++nl) {
137 for (nc=0; nc<ncols; ++nc) {
138 mask[nl][nc] = 1;
139 }
140 }
141
142
143 printf("============ Euclidean distance matrix between genes ============\n");
144 for (i = 0; i < ncols; i++) weight[i] = 1.0;
145 distMatrix = distancematrix(nrows, ncols, data, mask, weight, 'e', 0); // ZIAD: SIGKILL
146 if (!distMatrix)
147 { printf ("Insufficient memory to store the distance matrix\n");
148 free(weight);
149 return NULL;
150 }
151 /*printf(" Gene:");
152 for(i=0; i<nrows-1; i++) printf("%6d", i);
153 printf("\n");
154 for(i=0; i<nrows; i++)
155 { printf("Gene %2d:",i);
156 for(j=0; j<i; j++) printf(" %5.2f",distMatrix[i][j]);
157 printf("\n");
158 }
159 printf("\n");
160 */
161
162
163 for (ii=0;ii<nrows;++ii) {
164 if (mask[ii]) free(mask[ii]);
165 }
166 free(mask);
167 free(weight);
168 return distMatrix;
169 }
170
171
172 /* ========================================================================= */
173
example_hierarchical(int nrows,int ncols,double ** data,char * jobname,int k,double ** distmatrix)174 void example_hierarchical(int nrows, int ncols, double** data, char* jobname, int k, double** distmatrix)
175 /* Perform hierarchical clustering ... , double** distmatrix */
176 { int i, ii, nl, nc;
177 const int nnodes = nrows-1;
178 double* weight = malloc(ncols*sizeof(double));
179 int* clusterid;
180 Node* tree;
181 int** mask = NULL;
182 char* filename;
183 //char* filename2;
184
185
186 mask = (int **)calloc(sizeof(int*), nrows);
187 for (ii=0;ii<nrows;++ii) {
188 mask[ii] = (int *)calloc(sizeof(int),ncols);
189 }
190
191 for (nl=0; nl<nrows; ++nl) {
192 for (nc=0; nc<ncols; ++nc) {
193 mask[nl][nc] = 1;
194 }
195 }
196
197
198 for (i = 0; i < ncols; i++) weight[i] = 1.0;
199 printf("\n");
200
201
202 FILE *out1;
203
204 int n = 1 + strlen(jobname) + strlen("_C") + strlen(".ext");
205
206 if (k)
207 { int dummy = k;
208 do n++; while (dummy/=10);
209 }
210
211
212 filename = malloc(n*sizeof(char));
213
214 sprintf (filename, "%s_C%d.hie", jobname, k);
215 out1 = fopen( filename, "w" );
216
217 /*FILE *out2;
218 filename2 = malloc(n*sizeof(char));
219
220 sprintf (filename2, "%s_C%d.hi1", jobname, k);
221 out2 = fopen( filename2, "w" );*/
222
223 //HERE SHOULD USE method instead of 'xxx' (s,m,a,c)
224
225
226 printf("================ Pairwise single linkage clustering ============\n");
227 /* Since we have the distance matrix here, we may as well use it. */
228 tree = treecluster(nrows, ncols, 0, 0, 0, 0, 'e', 's', distmatrix);
229 /* The distance matrix was modified by treecluster, so we cannot use it any
230 * more. But we still need to deallocate it here.
231 * The first row of distmatrix is a single null pointer; no need to free it.
232 */
233 for (i = 1; i < nrows; i++) free(distmatrix[i]);
234 free(distmatrix);
235 if (!tree)
236 { /* Indication that the treecluster routine failed */
237 printf ("treecluster routine failed due to insufficient memory\n");
238 free(weight);
239 return;
240 }
241
242 #if 0
243 /* Andrej: This block looked like it was commented out
244 I took out some of the * / because they
245 were generating warning and blocked out the
246 entire section with #if 0 . */
247
248 /*fprintf(out2,"Node Item 1 Item 2 Distance\n");
249 for(i=0; i<nnodes; i++)
250 fprintf(out2,"%3d:%9d%9d %g\n",
251 -i-1, tree[i].left, tree[i].right, tree[i].distance);
252 printf("\n");
253 fclose(out2);*/
254 //free(tree);
255
256 /*
257 printf("================ Pairwise maximum linkage clustering ============\n");
258 tree = treecluster(nrows, ncols, data, mask, weight, 0, 'e', 'm', 0); */
259 /* Here, we let treecluster calculate the distance matrix for us. In that
260 * case, the treecluster routine may fail due to insufficient memory to store
261 * the distance matrix. For the small data sets in this example, that is
262 * unlikely to occur though. Let's check for it anyway:
263 */
264 if (!tree)
265 { /* Indication that the treecluster routine failed */
266 printf ("treecluster routine failed due to insufficient memory\n");
267 free(weight);
268 return;
269 }
270 printf("Node Item 1 Item 2 Distance\n");
271 for(i=0; i<nnodes; i++)
272 printf("%3d:%9d%9d %g\n",
273 -i-1, tree[i].left, tree[i].right, tree[i].distance);
274 printf("\n");
275 free(tree);
276
277
278
279 printf("================ Pairwise average linkage clustering ============\n");
280 tree = treecluster(nrows, ncols, data, mask, weight, 0, 'e', 'a', 0);
281 if (!tree)
282 { /* Indication that the treecluster routine failed */
283 printf ("treecluster routine failed due to insufficient memory\n");
284 free(weight);
285 return;
286 }
287 printf("Node Item 1 Item 2 Distance\n");
288 for(i=0; i<nnodes; i++)
289 printf("%3d:%9d%9d %g\n",
290 -i-1, tree[i].left, tree[i].right, tree[i].distance);
291 printf("\n");
292 free(tree);
293
294
295
296 printf("================ Pairwise centroid linkage clustering ===========\n");
297 tree = treecluster(nrows, ncols, data, mask, weight, 0, 'e', 'c', 0);
298 if (!tree)
299 { /* Indication that the treecluster routine failed */
300 printf ("treecluster routine failed due to insufficient memory\n");
301 free(weight);
302 return;
303 }
304 printf("Node Item 1 Item 2 Distance\n");
305 for(i=0; i<nnodes; i++)
306 printf("%3d:%9d%9d %g\n",
307 -i-1, tree[i].left, tree[i].right, tree[i].distance);
308 printf("\n");
309
310
311 #endif
312
313
314 printf("=============== Cutting a hierarchical clustering tree ==========\n");
315 clusterid = malloc(nrows*sizeof(int));
316 printf(" number of clusters %d \n",k);
317 cuttree (nrows, tree, k, clusterid);
318 for(i=0; i<nrows; i++)
319 fprintf(out1, "%09d\t%2d\n", i, clusterid[i]);
320 fprintf(out1, "\n");
321 fclose(out1);
322
323
324 for (ii=0;ii<nrows;++ii) {
325 if (mask[ii]) free(mask[ii]);
326 }
327 free(mask);
328 free(tree);
329 free(clusterid);
330 free(weight);
331 return;
332 }
333
getvoxlclusterdist(int * count,double ** cdata,int * clusterid,double ** data,char * jobname,int nclusters,int nrows,int ncols)334 void getvoxlclusterdist(int* count, double** cdata,
335 int* clusterid, double** data, char* jobname,
336 int nclusters, int nrows, int ncols)
337 {
338 int i, j, n;
339 char* filename4;
340 FILE *out4=NULL;
341 double* vcdata = malloc(nclusters*sizeof(double*));
342 double difference, difference1;
343
344 n = 1 + strlen(jobname) + strlen("_K_G") + strlen(".ext");
345
346 int dummy = nclusters;
347 do n++; while (dummy/=10);
348
349
350 filename4 = malloc(n*sizeof(char));
351 sprintf (filename4, "%s_K_G%d.vcd", jobname, nclusters);
352 out4 = fopen( filename4, "w" );
353
354 for (i = 0; i < nrows; i++){
355 difference = 0;
356 difference1 = 0;
357 for (j = 0; j < ncols; j++) {
358 difference1 = cdata[clusterid[i]][j]-data[i][j];
359 difference = difference + difference1*difference1;
360 }
361 vcdata[i] = sqrt(difference);
362 }
363
364 printf ("------- writing Cluster assignments + distances from centroids to file:\t\t"
365 " %s_K_G%d.vcd\n",jobname, nclusters);
366 for (i = 0; i < nrows; i++)
367 fprintf (out4, "%09d\t%d\t%7.3f\n", i, clusterid[i], vcdata[i]);
368 fclose(out4); out4=NULL;
369
370 /*for (i = 0; i < nrows; i++){
371 free(vcdata[i]);
372 }*/
373 free(vcdata);
374 return;
375
376 }
377
378
379
380 /* ========================================================================= */
381
382
example_kmeans(int nrows,int ncols,double ** data,int nclusters,int npass,char dist,char * jobname)383 void example_kmeans( int nrows, int ncols,
384 double** data,
385 int nclusters, int npass,
386 char dist, char* jobname)
387
388 /* Perform k-means clustering on genes */
389 {
390 int i, j, ii, nl, nc;
391 //const int nclusters = 3;
392 const int transpose = 0;
393 //const char dist = 'e';
394 const char method = 'a';
395 /* For method=='a', the centroid is defined as the mean over all elements
396 belonging to a cluster for each dimension.
397 For method=='m', the centroid is defined as the median over all elements
398 belonging to a cluster for each dimension.
399 */
400
401 //int npass = 1;
402 int ifound = 0;
403 int test=0;
404 double error;
405 double distance;
406 int** index;
407 int* count;
408 double* weight = malloc(ncols*sizeof(double));
409 int* clusterid = malloc(nrows*sizeof(int));
410 double** cdata = malloc(nclusters*sizeof(double*));
411 int** cmask = malloc(nclusters*sizeof(int*));
412 int** mask = NULL;
413 int n=0;
414 char* filename;
415 char* filename2;
416 char* filename3;
417 FILE *out1=NULL;
418 FILE *out2=NULL;
419 FILE *out3=NULL;
420
421 for (i = 0; i < nclusters; i++)
422 { cdata[i] = malloc(ncols*sizeof(double));
423 cmask[i] = malloc(ncols*sizeof(int));
424 }
425 for (i = 0; i < ncols; i++) weight[i] = 1.0;
426
427 mask = (int **)calloc(sizeof(int*), nrows);
428 for (ii=0;ii<nrows;++ii) {
429 mask[ii] = (int *)calloc(sizeof(int),ncols);
430 }
431
432 for (nl=0; nl<nrows; ++nl) {
433 for (nc=0; nc<ncols; ++nc) {
434 mask[nl][nc] = 1;
435 }
436 }
437
438
439 n = 1 + strlen(jobname) + strlen("_K_G") + strlen(".ext");
440
441 if (dist)
442 { int dummy = nclusters;
443 do n++; while (dummy/=10);
444 }
445
446 //avovk
447 printf("a je u omari :) \n");
448 filename = malloc(n*sizeof(char));
449 filename2 = malloc(n*sizeof(char));
450 filename3 = malloc(n*sizeof(char));
451 sprintf (filename, "%s_K_G%d.kgg", jobname, nclusters);
452 out1 = fopen( filename, "w" );
453 sprintf (filename2, "%s_K_G%d.dis", jobname, nclusters);
454 out2 = fopen( filename2, "w" );
455 sprintf (filename3, "%s_K_G%d.cen", jobname, nclusters);
456 out3 = fopen( filename3, "w" );
457
458 printf("======================== k-means clustering"
459 " ========================\n");
460
461 printf ("\n");
462 printf ("----- doing %d passes... go stretch your legs...\n",npass);
463 //npass = 3;
464 kcluster(nclusters,nrows,ncols,data,mask,weight,transpose,npass,method,dist,
465 clusterid, &error, &ifound);
466 printf ("Solution found %d times; ", ifound);
467 printf ("within-cluster sum of distances is %f\n", error);
468 printf ("------- writing Cluster assignments to file:\t\t"
469 " %s_K_G%d.kgg\n",jobname, nclusters);
470 for (i = 0; i < nrows; i++)
471 fprintf (out1, "%09d\t %d\n", i, clusterid[i]);
472 fclose(out1); out1=NULL;
473
474 printf ("------- writing Distance between clusters to file:\t %s_K_G%d.dis \n", jobname, nclusters);
475 fprintf (out2,"------- Distance between clusters:\n");
476 index = malloc(nclusters*sizeof(int*));
477 count = malloc(nclusters*sizeof(int));
478 for (i = 0; i < nclusters; i++) count[i] = 0;
479 for (i = 0; i < nrows; i++) count[clusterid[i]]++;
480 for (i = 0; i < nclusters; i++) index[i] = malloc(count[i]*sizeof(int));
481 for (i = 0; i < nclusters; i++) count[i] = 0;
482 for (i = 0; i < nrows; i++)
483 { int id = clusterid[i];
484 index[id][count[id]] = i;
485 count[id]++;
486 }
487
488 for (i = 0; i < nclusters-1; i++)
489 {
490 for (j = 1+i; j < nclusters; j++)
491 {
492 distance = clusterdistance(nrows, ncols, data, mask, weight, count[i], count[j], index[i], index[j], 'e', 'a', 0);
493 fprintf(out2,"Distance between %d and %d: %7.3f\n", i, j, distance);
494 // fprintf(stderr,"Distance between %d and %d: %7.3f\n", i, j, distance);
495 }
496 }
497
498 fclose(out2); out2=NULL;
499
500
501
502 printf ("------- writing Cluster centroids to file:\t\t%s_K_G%d.cen\n",jobname, nclusters);
503 fprintf (out3,"------- Cluster centroids:\n");
504 getclustercentroids(nclusters, nrows, ncols, data, mask, clusterid, cdata, cmask, 0, 'a');
505 fprintf(out3," coefficients:");
506 for(i=0; i<ncols; i++) fprintf(out3,"\t%7d", i);
507 fprintf(out3,"\n");
508 for (i = 0; i < nclusters; i++){
509 fprintf(out3,"Cluster %2d:", i);
510 for (j = 0; j < ncols; j++) fprintf(out3,"\t%7.3f", cdata[i][j]);
511 fprintf(out3,"\n");
512 }
513 fclose(out3); out3=NULL;
514 printf("Done...\n");
515
516 /* call function to calculate distance between each voxel and centroid */
517 /* we will need:
518 count - number of elements in cluster as we allready have it
519 cdata - cluster centroids
520 clusterid
521 data */
522
523
524 getvoxlclusterdist(count, cdata, clusterid, data, jobname,
525 nclusters, nrows, ncols);
526
527
528 for (i = 0; i < nclusters; i++) free(index[i]);
529 free(index);
530 free(count);
531
532 for (i = 0; i < nclusters; i++){
533 free(cdata[i]);
534 free(cmask[i]);
535 }
536 for (ii=0;ii<nrows;++ii) {
537 if (mask[ii]) free(mask[ii]);
538 }
539
540 free(cdata);
541 free(cmask);
542 free(clusterid);
543 free(weight);
544 free(mask);
545 return;
546 }
547
548
549 /* ========================================================================= */
550
main(int argc,char ** argv)551 int main(int argc, char **argv)
552 {
553 int ii=0, ncol=0, nrow=0, nl=0, nc=0, posi=0, posj=0, posk=0;
554 //int nclust=atoi(argv[2]);
555 MRI_IMAGE *im = NULL;
556 double *dar = NULL;
557 double **D = NULL;
558 //int **mask = NULL;
559 double** distmatrix = NULL;
560 //from command.c
561
562 int i = 1;
563 char* filename = 0;
564 char* jobname = 0;
565 int l = 0;
566 int k = 0;
567 int kh = 3;
568 int r = 1;
569 int s = 0;
570 int x = 2;
571 int y = 1;
572 int Rows, Columns;
573 char distmetric = 'u';
574 char arraymetric = '\0';
575 char method = 'm';
576 char cg = '\0';
577 char ca = '\0';
578 int ng = 0;
579 int na = 0;
580 while (i < argc)
581 { const char* const argument = argv[i];
582
583
584 i++;
585 if (strlen(argument)<2)
586 { printf("ERROR: missing argument\n");
587 return 0;
588 }
589 if (argument[0]!='-')
590 { printf("ERROR: unknown argument\n");
591 return 0;
592 }
593 if(!strcmp(argument,"--version") || !strcmp(argument,"-v"))
594 { display_version();
595 return 0;
596 }
597 if(!strcmp(argument,"--help") || !strcmp(argument,"-h"))
598 { display_help();
599 return 0;
600 }
601 if(!strcmp(argument,"-cg"))
602 { if (i==argc || strlen(argv[i])>1 || !strchr("am",argv[i][0]))
603 { printf ("Error reading command line argument cg\n");
604 return 0;
605 }
606 cg = argv[i][0];
607 i++;
608 continue;
609 }
610 if(!strcmp(argument,"-ca"))
611 { if (i==argc || strlen(argv[i])>1 || !strchr("am",argv[i][0]))
612 { printf ("Error reading command line argument ca\n");
613 return 0;
614 }
615 ca = argv[i][0];
616 i++;
617 continue;
618 }
619 if(!strcmp(argument,"-ng"))
620 { ng = 1;
621 continue;
622 }
623 if(!strcmp(argument,"-na"))
624 { na = 1;
625 continue;
626 }
627 switch (argument[1])
628 { case 'l': l=1; break;
629 case 'u':
630 { if (i==argc)
631 { printf ("Error reading command line argument u: "
632 "no job name specified\n");
633 return 0;
634 }
635 jobname = setjobname(argv[i],0);
636 i++;
637 break;
638 }
639 case 'f':
640 { if (i==argc)
641 { printf ("Error reading command line argument f: "
642 "no file name specified\n");
643 return 0;
644 }
645 filename = argv[i];
646 i++;
647 break;
648 }
649 case 'g':
650 { int g;
651 if (i==argc)
652 { printf ("Error reading command line argument g: parameter missing\n");
653 return 0;
654 }
655 g = readnumber(argv[i]);
656 if (g < 0 || g > 9)
657 { printf ("Error reading command line argument g: "
658 "should be between 0 and 9 inclusive\n");
659 return 0;
660 }
661 i++;
662 distmetric = getmetric(g);
663 break;
664 }
665
666 case 'k':
667 { if (i==argc)
668 { printf ("Error reading command line argument k: "
669 "parameter missing\n");
670 return 0;
671 }
672 k = readnumber(argv[i]);
673 if (k < 1)
674 { printf ("Error reading command line argument k: "
675 "a positive integer is required\n");
676 return 0;
677 }
678 i++;
679 break;
680 }
681
682 case 'c':
683 { if (i==argc)
684 { printf ("Error reading command line argument h: parameter missing\n");
685 return 0;
686 }
687 kh = readnumber(argv[i]);
688 if (kh < 1)
689 { printf ("Error reading command line argument h: a positive integer is required\n");
690 return 0;
691 }
692 i++;
693 break;
694 }
695 case 'r':
696 { if (i==argc)
697 { printf ("Error reading command line argument r: parameter missing\n");
698 return 0;
699 }
700 r = readnumber(argv[i]);
701 if (r < 1)
702 { printf ("Error reading command line argument r: "
703 "a positive integer is required\n");
704 return 0;
705 }
706 i++;
707 break;
708 }
709 case 'm':
710 { if (i==argc || strlen(argv[i])>1 || !strchr("msca",argv[i][0]))
711 { printf ("Error reading command line argument m: should be 'm', 's', 'c', or 'a'\n");
712 return 0;
713 }
714 method = argv[i][0];
715 i++;
716 break;
717 }
718
719
720 default: printf ("Unknown option\n");
721 }
722 }
723
724 if(jobname==0) jobname = setjobname(filename,1);
725
726 /* else
727 { display_help();
728 return 0;
729 }
730 */
731
732 if (argc < 2) {
733 display_help();
734 return 0;
735 }
736
737 //printf("num of clusters %d \n",nclust);
738 fprintf(stderr,"Patience, reading %s...\n ", filename);
739 im = mri_read_double_1D (filename);
740 /* ZIAD I get this warning
741 Aclustering.c:408: warning: passing argument 1 of mri_read_double_1D
742 discards qualifiers from pointer target type
743 Andrej: filename was declared as (const char *),
744 but the function expects (char *) */
745 if (!im) {
746 fprintf(stderr,"Error: Failed to read matrix data from %s\n",
747 argv[2]);
748 return(-1);
749 }
750 ncol = im->ny;
751 nrow = im->nx;
752 fprintf (stderr,"Have %d cols, %d rows\nNow...\n", ncol, nrow);
753
754 /* now just get the array and kill the rest */
755 dar = MRI_DOUBLE_PTR(im);
756
757 /* make sure that pointer is set to NULL in im, or risk hell */
758 mri_clear_data_pointer(im);
759 if (im) mri_free(im); im = NULL; /* now kill im */
760
761
762 /* for double loop*/
763 D = (double **)calloc(sizeof(double*), nrow-1);
764
765 for (ii=0;ii<(nrow-1);++ii) {
766 D[ii] = (double *)calloc(sizeof(double), ncol-1);
767 }
768
769 for (nl=1; nl<nrow; ++nl) {
770 for (nc=1; nc<ncol; ++nc) {
771 D[nl-1][nc-1] = dar[nl+nc*nrow];
772 }
773 //fprintf(stdout,"\n");
774 }
775
776 //show_data(nrows, ncols, data, mask);
777 //example_mean_median(nrows, ncols, data, mask);
778 //distmatrix = example_distance_gene(nrows, ncols, data, mask);
779 //if (distmatrix) example_hierarchical(nrows, ncols, data, mask, distmatrix);
780 //example_distance_array(nrows, ncols, data, mask);
781 if(k>0) example_kmeans(nrow-1, ncol-1, D, k, r, distmetric, jobname);
782
783 else
784 {
785 distmatrix = example_distance_gene(nrow-1, ncol-1, D); // ZIAD: goes2 SIGKILL error
786 if (distmatrix)
787 example_hierarchical(nrow-1, ncol-1, D, jobname, kh, distmatrix);
788 }
789
790
791
792 free(dar); dar = NULL; /* done with input array */
793 // To free D
794 for (ii=0;ii<(nrow-1);++ii) {
795
796 if (D[ii]) free(D[ii]);
797 }
798 free(D);
799 free(jobname);
800 //free();
801
802 fprintf (stderr,"\n");
803 return 0;
804 }
805
806 /* How to compile:
807 I copied Aclustering.c and cluster.c to /usr/local/bin/afni/src
808 then compile with:
809 gcc -Wall -Wno-unused Aclustering2.c cluster.c -o Aclustering -Inifti/niftilib -Inifti/nifticdf -Inifti/znzlib -L. -bind_at_load -l3DEdge -lmri -lf2c -lmri /usr/lib64/libXm.a -lXm -lXmu -lXp -lXpm -lXext -lXt -lX11 -lz -lexpat -lm -lc
810
811 The output are 3 files
812 - cluster_out.kgg -- as before, index number and cluster number
813 - cluster_out.dis --
814 - cluster_out.cen --
815
816
817 actually
818 cp /home/drejc/c_temp/Aclustering.c ./
819
820 gcc -Wall -g -Wno-unused Aclustering.c cluster.c -o Aclustering -Inifti/niftilib -Inifti/nifticdf -Inifti/znzlib -L. -l3DEdge -lmri -lf2c -lmri /usr/lib64/libXm.a -lXm -lXmu -lXp -lXpm -lXext -lXt -lX11 -lz -lexpat -lm -lc
821
822 cp Aclustering /home/drejc/segm4ziad/toandrej/out_mprage_10_smpl/
823
824 */
825