#ifndef lint
static char sccsid[]="@(#)spm_bias_mex.c 2.5 John Ashburner 03/05/12";
#endif
#include "spm_mapping.h"
#include "spm_datatypes.h"
/* #include "spm_sys_deps.h" */
/* #include "spm_vol_utils.h" */
#include
#include "mex.h"
#define EPS 2.2204460492503130808e-16
#define MAXB 128 /* Maximum number of bases functions in each dimension */
#define IGNORE_ZEROS
/********************************************************************************/
/* Extract a slice from a mapped volume
* vol - mapped volume
* i - slice number (first slice is 0)
* dat - extracted slice
*
* required: vol, i
* modified: dat
*/
static void get_slice(MAPTYPE *vol, int i, double dat[])
{
static double mat[16]={1,0,0,0, 0,1,0,0, 0,0,1,0, 0,0,0,1};
mat[14] = i+1.0;
slice(mat,dat,vol->dim[0],vol->dim[1], vol, 0,0);
if ((vol->dtype == SPM_UNSIGNED_CHAR) || (vol->dtype == SPM_SIGNED_SHORT) ||
(vol->dtype == SPM_SIGNED_INT) || (vol->dtype == SPM_SIGNED_SHORT_S) ||
(vol->dtype == SPM_SIGNED_INT_S) || (vol->dtype == SPM_SIGNED_CHAR) ||
(vol->dtype == SPM_UNSIGNED_SHORT) || (vol->dtype == SPM_UNSIGNED_INT) ||
(vol->dtype == SPM_UNSIGNED_SHORT_S) || (vol->dtype == SPM_UNSIGNED_INT_S))
{
/* Only a relatively small finite number of values can be stored as
integers. This can cause aliasing problems in the histograms.
The ideal solution is to assume an intensity PDF for each voxel
that is a top-hat, where the width is the scalefactor.
Adding to the histogram should really use the linear interpolation
kernel convolved with a top-hat of the same width as the scalefactor,
Instead, I've added some pseudo-random noise - much easier to implement
slightly faster, but not as good */
static float ran[] = {
0.141571,-0.193651,0.160932,-0.141983,0.438199,-0.012333,-0.40901,0.173834,
0.0148803,-0.278421,0.225009,-0.431753,0.464124,-0.292343,-0.338882,0.138221,
-0.499772,-0.164367,-0.2249,-0.455472,-0.406104,-0.0900026,0.316892,0.370517,
-0.477445,0.227177,0.348009,0.2286,0.455099,0.156351,0.242305,-0.155034,
0.384022,-0.152756,-0.440523,0.218415,0.458214,-0.343165,-0.0836469,-0.405965,
-0.0500535,0.369152,-0.108383,-0.247216,-0.145618,0.242978,0.150832,0.439793,
0.332799,-0.0300221,0.129866,-0.441812,0.0421874,-0.044274,0.363087,0.355197,
-0.0277443,0.286924,0.155982,-0.49996,-0.368763,-0.00512524,-0.461667,-0.272564,
-0.172117,0.399469,-0.18627,-0.248324,-0.0670109,0.342382,-0.315511,0.00817921,
-0.0477603,-0.174416,-0.119924,0.38648,0.261261,0.383766,-0.0425937,0.299202,
-0.365923,-0.434686,-0.124855,-0.126477,-0.0159776,0.469459,-0.157939,-0.247311,
0.0848869,0.0237036,-0.336581,-0.0136019,-0.00393925,0.343194,0.306198,0.357786,
0.109754,0.0657304,0.111899,-0.397023,-0.341684,-0.0863501,0.0604104,-0.231323,
0.284254,-0.112129,-0.469016,0.0855018,0.0585585,-0.299304,-0.412578,0.43323,
-0.24062,-0.295829,-0.450792,0.106161,0.0463487,-0.404163,0.136996,-0.0570517,
-0.433618,-0.125707,-0.250897,0.424875,0.129499,0.378309,0.141674,0.298391,
-0.064974,0.48114,-0.404042,0.0274824,0.0456458,-0.215657,-0.129197,-0.435307,
0.0448091,0.336376,-0.354678,-0.32848,-0.431953,0.324012,-0.366029,0.384786,
0.0147374,0.463636,-0.379505,-0.45171,-0.119848,-0.0872088,-0.0986087,-0.0790029,
-0.123046,0.407337,0.170162,0.461839,-0.337021,0.248649,-0.125934,-0.0457635,
-0.461439,0.062432,-0.127688,0.292784,0.295231,-0.117086,-0.24721,-0.157072,
0.467804,-0.0201916,-0.131672,0.264567,-0.122851,0.400306,-0.316568,-0.131683,
0.417457,0.0159161,-0.409693,0.235311,-0.495288,0.103123,0.456867,-0.102568,
0.231551,0.184639,0.478503,-0.296215,0.0933028,0.451563,-0.239728,0.0146776,
0.136332,-0.0990409,-0.0133887,0.250458,-0.373803,-0.456938,-0.129061,0.193302,
0.435825,-0.0224244,-0.3709,-0.0161688,0.445597,-0.132256,-0.171514,0.272875,
-0.202726,-0.322146,0.190799,-0.236055,-0.0422609,0.343692,0.381499,0.20002,
0.255691,0.474515,-0.0978473,-0.368721,0.224734,0.399518,-0.329299,-0.456971,
-0.0208408,-0.406064,0.15005,0.452278,-0.0422873,0.0368806,-0.433513,-0.00612992,
-0.082459,-0.207743,-0.210336,0.253846,-0.403204,-0.423083,0.220915,0.264912,
0.157945,0.310409,-0.125757,-0.193769,-0.1293,0.206747,-0.331634,0.313721,
-0.0337717,0.222286,0.494869,-0.137497,0.230828,0.149667,0.181339,-0.492388,
0.154149,0.445235,0.113271,0.282928,-0.496847,0.296958,0.141817,-0.321521,
0.0294007,-0.281256,0.0480524,-0.441763,0.0875871,-0.0838964,-0.313555,-0.436089,
-0.425217,-0.189963,0.444085,0.480727,0.0551155,0.488526,0.19156,-0.258339,
0.309814,0.434512,-0.371197,0.186827,-0.202755,0.147221,-0.0362412,0.42279,
-0.258278,0.160157,-0.167713,-0.023247,-0.0312267,0.205929,-0.260051,0.217203,
0.365197,-0.0890103,-0.0752128,0.454286,0.381419,0.198588,-0.194599,0.329326,
0.470605,-0.199897,0.498062,-0.0592516,-0.493773,-0.208099,0.191878,-0.00715073,
-0.416522,-0.304213,0.477595,-0.133801,-0.360568,-0.485244,0.140641,0.237722,
-0.473243,-0.397847,0.39333,0.263876,0.0826645,0.18543,0.467134,0.0503444,
0.277197,0.11035,0.488593,-0.451738,0.485389,-0.295264,0.412477,0.165546,
-0.0376786,-0.451667,-0.0395608,0.299951,-0.210567,0.195102,-0.240721,0.21323,
0.22037,0.233277,0.122331,0.489779,-0.347622,-0.296677,0.319317,-0.441565,
0.0385068,-0.309847,0.0994808,-0.207742,-0.408712,0.00676728,0.384132,0.115568,
-0.453608,0.451924,-0.33097,0.326671,0.111377,0.347299,-0.385877,0.14922,
-0.385185,-0.0265878,0.183181,-0.366701,-0.0359218,-0.42875,0.0812184,0.0659789,
-0.244738,-0.261506,-0.484046,-0.115255,0.2573,0.07516,-0.0919145,-0.304292,
0.0121743,0.213346,0.367386,0.241805,0.494833,0.366748,-0.0141872,0.403305,
-0.479686,0.136787,0.444056,-0.45922,-0.214328,0.177254,0.0861783,-0.43503};
int j, iran = i*13;
double sf = vol->scale[i];
iran = i*10;
for(j=0;jdim[0]*vol->dim[1];j++)
{
iran = (iran+1)%397;
#ifdef IGNORE_ZEROS
if (dat[j])
#endif
dat[j] += ran[iran]*sf;
}
}
}
/********************************************************************************/
/* Normalise the histogram so that it integrates to unity, taking the width of
* the bins into account.
* nh - number of histogram bins
* h - histogram
* binwidth - width of each histogram bin (1/scale)
*
* required: nh, h, binwidth
* modified: h
*/
static void normhist(int nh, double h[], double binwidth, int *psh)
{
double sh = 0.0;
int i;
for(i=0; i0.0)
#else
if (y0>=0.0)
#endif
{
int iy;
double dy, sc = 0.0;
for(i0=0;i0=0 && y<(nh-1))
{
iy = floor(y);
dy = y-iy;
h[iy] += 1-dy;
h[iy+1] += dy;
}
}
}
}
}
normhist(nh,h,1/s,psh);
}
/********************************************************************************/
/********************************************************************************/
/* Create voxel values in weighting images and increment likelihoods using
* lookup tables
* dat - raw voxel intensity
* sc - estimated scaling
* nh - number of bins in lookup table
* s - scale from image intensity to bin number
* wt0 - first weight
* wt1 - second weight
* ll - incremented log-likelihood
*
* required: dat, sc, nh, s
* modified: wt0, wt1, ll
*/
static void weights(double dat, double sc,
int nh, double s, double h[],
double *wt0, double *wt1, double *ll)
{
double y, dy, f, df, t;
int iy;
y = dat*sc*s;
#ifdef IGNORE_ZEROS
if (y>0.0 && y=0.0 && ydim[0];
n1 = vol->dim[1];
n2 = vol->dim[2];
m0 = mxGetN(prhs[1]);
m1 = mxGetN(prhs[2]);
m2 = mxGetN(prhs[3]);
if (m0>MAXB || m1>MAXB || m2>MAXB)
{
free_maps(vol, 1);
mexErrMsgTxt("Too many basis functions.");
}
B0 = mxGetPr(prhs[1]);
B1 = mxGetPr(prhs[2]);
B2 = mxGetPr(prhs[3]);
if (mxGetM(prhs[1])!=n0 || mxGetM(prhs[2])!=n1 || mxGetM(prhs[3])!=n2)
{
free_maps(vol, 1);
mexErrMsgTxt("Inappropriate matrix dimensions.");
}
plhs[0] = mxCreateDoubleMatrix(m0*m1*m2,1, mxREAL);
plhs[1] = mxCreateDoubleMatrix(1,1, mxREAL);
mxGetPr(plhs[1])[0] = constr(m0,m1,m2,n0,n1,n2,mxGetPr(plhs[0]),vol,B0,B1,B2);
free_maps(vol, 1);
}
/********************************************************************************/
void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
int m0, m1, m2, n0, n1, n2, nh, sh;
MAPTYPE *vol, *get_maps();
double *B0, *B1, *B2, mx, *f, *Alpha, *Beta, *h;
if (nrhs == 4 && nlhs <=2)
{
mexFunction1(nlhs, plhs, nrhs, prhs);
return;
}
if (nrhs != 6 || nlhs > 5) mexErrMsgTxt("Inappropriate usage.");
vol = get_maps(prhs[0], &n0);
if (n0!=1)
{
free_maps(vol, n0);
mexErrMsgTxt("Inappropriate usage.");
}
for(n0=1; n0dim[0];
n1 = vol->dim[1];
n2 = vol->dim[2];
m0 = mxGetN(prhs[1]);
m1 = mxGetN(prhs[2]);
m2 = mxGetN(prhs[3]);
if (mxGetM(prhs[1])!=n0 || mxGetM(prhs[2])!=n1 || mxGetM(prhs[3])!=n2
|| mxGetM(prhs[4])*mxGetN(prhs[4]) != m0*m1*m2
|| mxGetM(prhs[5])*mxGetN(prhs[5])!=2)
{
free_maps(vol, 1);
mexErrMsgTxt("Inappropriate matrix dimensions.");
}
if (m0>MAXB || m1>MAXB || m2>MAXB)
{
free_maps(vol, 1);
mexErrMsgTxt("Too many basis functions.");
}
B0 = mxGetPr(prhs[1]);
B1 = mxGetPr(prhs[2]);
B2 = mxGetPr(prhs[3]);
f = mxGetPr(prhs[4]);
mx = mxGetPr(prhs[5])[0];
nh = floor(mxGetPr(prhs[5])[1]);
if (nh<1 || nh > 10240)
{
free_maps(vol, 1);
mexErrMsgTxt("Inappropriate number of bins in histogram.");
}
plhs[0] = mxCreateDoubleMatrix(m0*m1*m2,m0*m1*m2, mxREAL);
plhs[1] = mxCreateDoubleMatrix(m0*m1*m2,1, mxREAL);
plhs[2] = mxCreateDoubleMatrix(1,1, mxREAL);
plhs[3] = mxCreateDoubleMatrix(nh,1, mxREAL);
plhs[4] = mxCreateDoubleMatrix(1,1, mxREAL);
/* Initialised to zero */
Alpha = mxGetPr(plhs[0]);
Beta = mxGetPr(plhs[1]);
h = mxGetPr(plhs[3]);
mxGetPr(plhs[2])[0] = iteration(m0,m1,m2, n0,n1,n2, Alpha,Beta, vol, B0,B1,B2,f,mx, nh,h, &sh);
mxGetPr(plhs[4])[0] = (double)sh;
free_maps(vol, 1);
}