← Back to branch summary

~nipy-developers/nipy/fff2

« back to all changes in this revision

Viewing changes to matlab/fffMat_clusters.c

  • Committer: Bertrand THIRION
  • Date: 2008-04-03 17:29:55 UTC
  • mfrom: (13.1.5 fff2)
  • Revision ID: bt206016@is143015-20080403172955-w7v1pdjdriofvzzs
Merged Alexis'changes

Show diffs side-by-side

added added

removed removed

Lines of Context:
matlab/fffMat_clusters.c
 
1
#include "fffMat.h"
 
2
 
 
3
#include <fff_graphlib.h>
 
4
#include <fff_field.h>
 
5
 
 
6
#include <string.h>
 
7
 
 
8
/* 
 
9
 
 
10
SYNTAX: [size, label, depth] = fffMat_clusters( XYZ, T ); 
 
11
 
 
12
Warning: the input array XYZ must be in C-compatible long int format
 
13
(most often, 32 bits or 64 bits).
 
14
 
 
15
The returned arrays (size, label, depth) are all C-compatible long int
 
16
arrays.
 
17
 
 
18
*/ 
 
19
 
 
20
void mexFunction( int nlhs, mxArray *plhs[], 
 
21
                  int nrhs, const mxArray*prhs[] )
 
22
{ 
 
23
  mxClassID mxLong;
 
24
  gsl_matrix_long XYZ; 
 
25
  gsl_vector* T;
 
26
  fff_graph* G = NULL; 
 
27
  long *label, *bufL, *size; 
 
28
  gsl_vector_long depth; 
 
29
  int nvox, edges, clusters, maxima; 
 
30
  int i, j; 
 
31
 
 
32
  /* Check for proper number of arguments. */
 
33
  if(nrhs < 2) {
 
34
    mexErrMsgTxt("Two input arguments required.");
 
35
  } 
 
36
  else if(nlhs > 3) {
 
37
    mexErrMsgTxt("Too many output arguments");
 
38
  }
 
39
  
 
40
  /* Input */ 
 
41
  nvox = mxGetN( prhs[0] ); 
 
42
 
 
43
  /* Consistency check */ 
 
44
  if ( ( mxGetM( prhs[0] )!= 3 ) ||
 
45
       ( mxGetN( prhs[1] )!= nvox ) )
 
46
    mexErrMsgTxt("Inconsistent input arrays.");
 
47
  
 
48
  /* Check input type */ 
 
49
  if ( ! mxArray_is_long( prhs[0], &mxLong ) )
 
50
    mexErrMsgTxt("Input array XYZ is not in proper int format.");
 
51
 
 
52
  /* XYZ matrix */ 
 
53
  XYZ.size1 = nvox; 
 
54
  XYZ.size2 = 3; 
 
55
  XYZ.tda = 3; 
 
56
  XYZ.data = (long*) mxGetPr( prhs[0] ); 
 
57
  XYZ.owner = 0; 
 
58
  XYZ.block = NULL; 
 
59
  
 
60
  /* Field values */ 
 
61
  T = gsl_vector_from_mxArray( prhs[1] ); 
 
62
  
 
63
  /* Output */ 
 
64
  plhs[1] = mxCreateNumericMatrix( nvox, 1, mxLong, mxREAL );
 
65
  plhs[2] = mxCreateNumericMatrix( nvox, 1, mxLong, mxREAL );
 
66
 
 
67
  /* Make a graph out of the input coordinates */
 
68
  edges = fff_graph_grid( &G, &XYZ, 18 );
 
69
  
 
70
  /* Assign cluster label to each voxel */ 
 
71
  label = (long*) mxGetPr( plhs[1] );
 
72
  clusters = fff_graph_cc_label( label, G );
 
73
 
 
74
  /* Compute the cluster sizes */ 
 
75
  plhs[0] = mxCreateNumericMatrix( clusters, 1, mxLong, mxREAL );
 
76
  size = (long*) mxGetPr( plhs[0] );
 
77
  bufL = label; 
 
78
  for( j=0; j<nvox; j++, bufL++ ) {
 
79
    i = *bufL; 
 
80
    size[i] += 1.0; 
 
81
    *bufL = i+1; /* To return indices starting from 1... */ 
 
82
  }
 
83
 
 
84
  /* Compute depth for each voxel */
 
85
  depth.size = (size_t) nvox; 
 
86
  depth.stride = 1; 
 
87
  depth.data = (long*) mxGetPr( plhs[2] ); 
 
88
  depth.owner = 0; 
 
89
  depth.block = NULL; 
 
90
  maxima = fff_field_maxima( &depth, G, T ); 
 
91
 
 
92
  /* Free memory */ 
 
93
  fff_graph_delete( G ); 
 
94
  gsl_vector_free( T ); 
 
95
 
 
96
  return; 
 
97
 
 
98
}
 
99
 
 
100