A better way to create a response function for CSD #404

ChantalTax · 2014-07-31T18:36:49Z

This implementation includes the recursive calibration of the response function for CSD, Using an FA threshold is not a very robust method.It is dependent on the dataset (non-informed used subjectivity), and still depends on the diffusion tensor (FA and first eigenvector), which has low accuracy at high b-value. This function recursively calibrates the response function, for more information see Tax et al. NeuroImage 2014.

This implementation is slow only on one line (788) which calls real_sph_harm. I would appreciate any help to optimize this function. Thx!

arokem · 2014-07-31T18:57:45Z

dipy/reconst/csdeconv.py

+        dirs = dirs[single_peak_mask]
+
+        for num_vox in range(0, data.shape[0]):
+            rotmat = vec2vec_rotmat(dirs[num_vox, 0], np.array([0, 0, 1]))


It looks like maybe some of these computations (e.g. this line) could be moved outside of the num_it loop (above line 762)

Or I might be misunderstanding, not sure.

Hey @arokem, thanks a lot for your comments! I'm not sure if this is the case? I want to rotate the main peak direction of every voxel towards the z-axis here. So I actually need to compute the rotation matrix that does this for every voxel independently (this is rotmat). Then I rotate the gradients accordingly, and fit spherical harmonics. Maybe I overlook something, could you perhaps comment on that? Thank you!

Yep - gotcha - sorry about that.

On Thu, Aug 21, 2014 at 10:53 AM, ChantalTax notifications@github.com
wrote:

In dipy/reconst/csdeconv.py:

csd_peaks = peaks_from_model(model=csd_model,

data=data,

sphere=sphere,

relative_peak_threshold=peak_thr,

min_separation_angle=25,

parallel=parallel)

dirs = csd_peaks.peak_dirs

vals = csd_peaks.peak_values

single_peak_mask = (vals[:, 1] / vals[:, 0]) < peak_thr

data = data[single_peak_mask]

dirs = dirs[single_peak_mask]

for num_vox in range(0, data.shape[0]):

rotmat = vec2vec_rotmat(dirs[num_vox, 0], np.array([0, 0, 1]))

Hey @arokem https://github.com/arokem, thanks a lot for your comments!
I'm not sure if this is the case? I want to rotate the main peak direction
of every voxel towards the z-axis here. So I actually need to compute the
rotation matrix that does this for every voxel independently (this is
rotmat). Then I rotate the gradients accordingly, and fit spherical
harmonics. Maybe I overlook something, could you perhaps comment on that?
Thank you!

—
Reply to this email directly or view it on GitHub
https://github.com/nipy/dipy/pull/404/files#r16555307.

arokem · 2014-07-31T19:00:16Z

Hey @ChantalTax - looks overall great, but I haven't done much of a fine-grained review yet. Would you mind rebasing this on top of master before we move on? Let me know if you need help with that

MrBago · 2014-07-31T19:15:33Z

dipy/reconst/csdeconv.py

+            x, y, z = rot_gradients[where_dwi].T
+            r, theta, phi = cart2sphere(x, y, z)
+            # for the gradient sphere
+            B_dwi = real_sph_harm(m, n, theta[:, None], phi[:, None])


it seems like B_dwi can be computed outside of the loop too,

Hey @MrBago, many thanks for your comments! As I wrote in the above comment, I am not sure if this is possible, but maybe I overlook something? I want to rotate the main peak direction of every voxel towards the z-axis here. So I actually need to compute the rotation matrix that does this for every voxel independently (this is rotmat). Then I rotate the gradients accordingly, and fit spherical harmonics.

Sorry @ChantalTax, I started shooting my mouth off before I gave this a proper look over. This looks really awesome overall. I think you're right, there is no way to avoid calling real_sh_harm inside of this loop, which is too bad because the scipy implementation is quite slow (but I hear that might be improving soon).

Looking over this, the main thing I notice is that you're not enforcing axial symmetry on the response function. I think that might be a good idea because you're using a vector alignment to get your rotation matrix, this leaves the axial rotation angle ambiguous (which makes sense for this application), but then I think you should force the response function to be axially symmetric (m != 0 harmonics have zero coefficient). This might also allow you to improve performance because you don't need to build the full B_dwi matrix every time. That being said, you've spent a lot more time on this than I have, so just let me know if I'm rambling nonsense.

In terms of performance, can i suggest using a HemiSphere instead of a full sphere. The HemiSphere implementation in dipy takes advantages of the antipodal symmetry inherent in diffusion imaging. It works everywhere a sphere is needed (but doesn't look as good for displaying ODFs/FODs and such), but is twice as fast. You can import the HemiSphere version of 'symmetric724' by using from dipy.reconst.peaks import default_sphere, or you can build a HemiSphere from any sphere by using hemisphere = HemiSphere.from_sphere(sphere). This will eliminate the colinear lines in your B_dwi matrix, which aren't helping you anyways.

Thanks for you great contribution, let me know if you need any help with my suggestions above.

Please see my suggestion above for a significantly faster way to
evaluate sph_harm, which may help out here.

@MrBago yes you are totally right about the axial symmetric response function. Up to now I didn't do this because, if I'm not mistaken, the function sh_to_rh does this before actual CSD. But it would be nicer to output an axially symmetric response function, I totally agree.

I'm not entirely sure what you mean by using a hemisphere, and to which lines of codes you refer to? Could you shed some light on that?
Thanks!

Garyfallidis · 2014-08-13T23:19:29Z

@stefanv I remember that you had identified in the past that real_sph_harm is too slow and you had some ideas on how to better implement this. For @ChantalTax the performance of real_sph_harm is crucial as she is calling it many many times. Would you mind suggesting a performance boost for this function. Make it or show to Chantal how she could do it?

stefanv · 2014-08-13T23:39:53Z

This will be fixed after this year's scipy GSoC.

Garyfallidis · 2014-08-14T13:13:26Z

Oh that is great @stefanv. Until we start depending on the latest scipy version we will need to have something working nicely in dipy. Would it be possible to let us know when that update is ready so that we can copy it temporarily in Dipy? Would that be reasonable or this will be C code?

stefanv · 2014-08-14T16:19:21Z

It may require access to Fortran. @pv, can you confirm the requirements
for the new, faster spherical harmonics calculation?

ken-sakaie · 2014-08-14T16:51:07Z

Eleftherios,

Couldn’t help but chime in on this one.

Numerical recipes has some really simple code for calculating these. See the chapter on special functions, sub-chapter “Spherical Harmonics”

The code can be stripped down even more to speed things up.

I’ve attached some code I use (in C) to calculate the real spherical harmonics.

Feel free to use it if it helps.

Sorry I don’t have time to do much coding.

Ken

From: Eleftherios Garyfallidis [mailto:notifications@github.com]
Sent: Wednesday, August 13, 2014 7:20 PM
To: nipy/dipy
Subject: Re: [dipy] A better way to create a response function for CSD (#404)

@stefanvhttps://github.com/stefanv I remember that you had identified in the past that real_sph_harm is too slow and you had some ideas on how to better implement this. For @ChantalTaxhttps://github.com/ChantalTax the performance of real_sph_harm is crucial as she is calling it many many times. Would you mind suggesting a performance boost for this function. Make it or show to Chantal how she could do it?

—
Reply to this email directly or view it on GitHubhttps://github.com//pull/404#issuecomment-52125000.

Please consider the environment before printing this e-mail

Cleveland Clinic is ranked as one of the top hospitals in America by U.S.News & World Report (2014).
Visit us online at http://www.clevelandclinic.org for a complete listing of our services, staff and locations.

Confidentiality Note: This message is intended for use only by the individual or entity to which it is addressed and may contain information that is privileged, confidential, and exempt from disclosure under applicable law. If the reader of this message is not the intended recipient or the employee or agent responsible for delivering the message to the intended recipient, you are hereby notified that any dissemination, distribution or copying of this communication is strictly prohibited. If you have received this communication in error, please contact the sender immediately and destroy the material in its entirety, whether electronic or hard copy.

Thank you.

/******************************************************************************
*

```
    FILE: ylm_calc.c
```

   USAGE: void ylm_calc(const int lmax, float z, float phi, float \* ylm);

 PURPOSE: Calculates real parts of spherical harmonic with antipodal

          symmetry: Ylm(z=cos(theta),phi)

   INPUT: 
    lmax: maximum degree, l. MUST BE EVEN.

      z: cosine of angle with respect to z axis. ASSUMES -1<=z<=1

```
      phi: azimuthal angle, IN RADIANS
```
OUTPUT:
ylm: array of spherical harmonics.

           has (lmax+1)*(lmax/2+1) components.

```
  Functions are in order:
```
```
  N=1;
```
```
  for L=0:LMAX
```
```
      ylm(N) = YLM(L,0);
```
```
      N=N+1;
```
```
      for M=1:L
```
```
          ylm(N) = -2*Im(YLM(L,M));
```
```
          N=N+1;
```
```
          ylm(N) = 2*Re(YLM(L,M));
```
```
          N=N+1;
```
```
      end
```
```
  end
```
Copyright 2009 by Ken Sakaie. Please don't sue me.
$URL: svn+ssh://pub1/home/sakaiek/repos/trunk/workspace/chardi/fod/ylm_calc.c $
$Author: sakaiek $
$Date: 2009-02-12 13:02:06 -0500 (Thu, 12 Feb 2009) $
$Rev: 1081 $
$Id: ylm_calc.c 1081 2009-02-12 18:02:06Z sakaiek $

*****************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define M_4PI 12.56637061435917295405

void ylm_calc( int lmax, float z, float phi, float * ylm)
{
/* allocate workspace /
float * work = (float)malloc((lmax+1)*sizeof(float));
if(work == NULL) {
fprintf(stderr, "ylm_calc(work): memory error\n");
exit(1);
}

/* coefficient for legendre functions */
int l;
for(l=0;l<=lmax;l++) work[l]=sqrt((2*l+1)/M_4PI);

/* variables */
int m;
float Pmm=1.0, mfac=1.0, lmfac=0.0;         /* for legendre */
float P_1=0.0, P_2=0.0, Pm1m=0.0, Plm=0.0;  /* for legendre recursion */
int j0=2,j=1;                   /* indices for output */
int meven=1, leven=1;               /* flags: check m,l even */
float aylm=0.0, sm=0.0, cm=0.0;     /* function values */

/* calculate */
for(m=0;m<=lmax;m++) {

    /* Ymm */
    if ( m == 0 ) {
        /* legendre */
        ylm[0]=work[0];
        /* set indices */
        meven=0;
    } else {
        /* phase factors */
        sm = -2.0 * sin( m * phi );
        cm =  2.0 * cos( m * phi );
        /* legendre */
        Pmm = -(2*m-1)*sqrt(1-z*z)*Pmm;
        mfac = mfac / sqrt( (m+m)*(m+m-1) );
        aylm = work[m] * mfac * Pmm;
        /* index, assign */
        if ( meven==1 ) {
            j = j0;
            ylm[j]   = aylm * sm;
            ylm[j+1] = aylm * cm;
            j0 += 2*m + 3;
            j  += 2*m + 1;
            meven=0;
        } else {
            meven=1;
        }
    }

    /* Ym+1,m */
    if ( m < lmax ) {
        /* legendre */
        Pm1m = z * (2*m+1) * Pmm;
        lmfac = mfac / sqrt(m+m+1);
        aylm = work[m+1] * lmfac * Pm1m;
        /* index */
        if ( meven==1 ) { /* m odd, m+1 even */
            j=j0;
            ylm[j]   = aylm * sm;
            ylm[j+1] = aylm * cm;
            j0 += 2;
            j  += 2*(m+1)+1;
        }
    }

    if ( m < (lmax-1) ) {
        P_2 = Pmm;
        P_1 = Pm1m;
        leven = 1 - meven;
        for(l=m+2;l<=lmax;l++) {
            /* factorials */
            lmfac *= sqrt(l-m) / sqrt(l+m);
            /* legendre */
            Plm = z*( 2*l-1 )*P_1 - (l+m-1)*P_2;
            Plm = Plm / (l-m);
            P_2 = P_1;
            P_1 = Plm;
            aylm = work[l] * lmfac * Plm;
            /* index */
            if ( leven == 1 ) {
                if ( m== 0 ) ylm[j] = aylm;
                else {
                    ylm[j]   = aylm * sm;
                    ylm[j+1] = aylm * cm;
                }
                j = j + 2*l + 1;
                leven=0;
            } else {
                leven=1;
            }
        }
    }
}

/* clean up */
free(work);

}
#undef M_4PI

stefanv · 2014-08-14T17:02:00Z

Just note that numerical recipes code is not licensed under a free license
exactly, although I imagine it is intended to be used so.

On Thu, Aug 14, 2014 at 6:51 PM, ken-sakaie notifications@github.com
wrote:

Eleftherios,

Couldn’t help but chime in on this one.

Numerical recipes has some really simple code for calculating these. See
the chapter on special functions, sub-chapter “Spherical Harmonics”

The code can be stripped down even more to speed things up.

I’ve attached some code I use (in C) to calculate the real spherical
harmonics.

Feel free to use it if it helps.

Sorry I don’t have time to do much coding.

Ken

From: Eleftherios Garyfallidis [mailto:notifications@github.com]
Sent: Wednesday, August 13, 2014 7:20 PM
To: nipy/dipy
Subject: Re: [dipy] A better way to create a response function for CSD
(#404)

@stefanvhttps://github.com/stefanv I remember that you had identified
in the past that real_sph_harm is too slow and you had some ideas on how to
better implement this. For @ChantalTaxhttps://github.com/ChantalTax the
performance of real_sph_harm is crucial as she is calling it many many
times. Would you mind suggesting a performance boost for this function.
Make it or show to Chantal how she could do it?

—
Reply to this email directly or view it on GitHub<
https://github.com/nipy/dipy/pull/404#issuecomment-52125000>.

Please consider the environment before printing this e-mail

Cleveland Clinic is ranked as one of the top hospitals in America by
U.S.News & World Report (2014).
Visit us online at http://www.clevelandclinic.org for a complete listing
of our services, staff and locations.

Confidentiality Note: This message is intended for use only by the
individual or entity to which it is addressed and may contain information
that is privileged, confidential, and exempt from disclosure under
applicable law. If the reader of this message is not the intended recipient
or the employee or agent responsible for delivering the message to the
intended recipient, you are hereby notified that any dissemination,
distribution or copying of this communication is strictly prohibited. If
you have received this communication in error, please contact the sender
immediately and destroy the material in its entirety, whether electronic or
hard copy.

Thank you.

/******************************************************************************
*

FILE: ylm_calc.c
*

USAGE: void ylm_calc(const int lmax, float z, float phi, float * ylm);
*

PURPOSE: Calculates real parts of spherical harmonic with antipodal

symmetry: Ylm(z=cos(theta),phi)
*

INPUT:
lmax: maximum degree, l. MUST BE EVEN.

z: cosine of angle with respect to z axis. ASSUMES -1<=z<=1

phi: azimuthal angle, IN RADIANS
*

OUTPUT:
ylm: array of spherical harmonics.

has (lmax+1)*(lmax/2+1) components.
*

Functions are in order:

N=1;

for L=0:LMAX

ylm(N) = YLM(L,0);

N=N+1;

for M=1:L

ylm(N) = -2*Im(YLM(L,M));

N=N+1;

ylm(N) = 2*Re(YLM(L,M));

N=N+1;

end

end
*
Copyright 2009 by Ken Sakaie. Please don't sue me.
$URL:
svn+ssh://pub1/home/sakaiek/repos/trunk/workspace/chardi/fod/ylm_calc.c $
$Author: sakaiek $
$Date: 2009-02-12 13:02:06 -0500 (Thu, 12 Feb 2009) $
$Rev: 1081 $
$Id: ylm_calc.c 1081 2009-02-12 18:02:06Z sakaiek $

*****************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define M_4PI 12.56637061435917295405

void ylm_calc( int lmax, float z, float phi, float * ylm)
{
/* allocate workspace /
float * work = (float)malloc((lmax+1)*sizeof(float));
if(work == NULL) {
fprintf(stderr, "ylm_calc(work): memory error\n");
exit(1);
}

/* coefficient for legendre functions _/
int l;
for(l=0;l<=lmax;l++) work[l]=sqrt((2_l+1)/M_4PI);

/* variables /
int m;
float Pmm=1.0, mfac=1.0, lmfac=0.0; / for legendre /
float P_1=0.0, P_2=0.0, Pm1m=0.0, Plm=0.0; / for legendre recursion /
int j0=2,j=1; / indices for output /
int meven=1, leven=1; / flags: check m,l even /
float aylm=0.0, sm=0.0, cm=0.0; / function values */

/* calculate */
for(m=0;m<=lmax;m++) {

/* Ymm /
if ( m == 0 ) {
/ legendre /
ylm[0]=work[0];
/ set indices /
meven=0;
} else {
/ phase factors /
sm = -2.0 * sin( m * phi );
cm = 2.0 * cos( m * phi );
/ legendre _/
Pmm = -(2_m-1)_sqrt(1-z_z)Pmm;
mfac = mfac / sqrt( (m+m)(m+m-1) );
aylm = work[m] * mfac * Pmm;
/* index, assign _/
if ( meven==1 ) {
j = j0;
ylm[j] = aylm * sm;
ylm[j+1] = aylm * cm;
j0 += 2_m + 3;
j += 2*m + 1;
meven=0;
} else {
meven=1;
}
}

/* Ym+1,m /
if ( m < lmax ) {
/ legendre _/
Pm1m = z * (2_m+1) * Pmm;
lmfac = mfac / sqrt(m+m+1);
aylm = work[m+1] * lmfac * Pm1m;
/* index /
if ( meven==1 ) { / m odd, m+1 even /
j=j0;
ylm[j] = aylm * sm;
ylm[j+1] = aylm * cm;
j0 += 2;
j += 2(m+1)+1;
}
}

if ( m < (lmax-1) ) {
P_2 = Pmm;
P_1 = Pm1m;
leven = 1 - meven;
for(l=m+2;l<=lmax;l++) {
/* factorials /
lmfac *= sqrt(l-m) / sqrt(l+m);
/ legendre /
Plm = z( 2_l-1 )_P_1 - (l+m-1)P_2;
Plm = Plm / (l-m);
P_2 = P_1;
P_1 = Plm;
aylm = work[l] * lmfac * Plm;
/ index _/
if ( leven == 1 ) {
if ( m== 0 ) ylm[j] = aylm;
else {
ylm[j] = aylm * sm;
ylm[j+1] = aylm * cm;
}
j = j + 2_l + 1;
leven=0;
} else {
leven=1;
}
}
}
}

/* clean up */
free(work);

}
#undef M_4PI

—
Reply to this email directly or view it on GitHub
#404 (comment).

pv · 2014-08-14T17:09:18Z

The revised sph_harm implementation is faster by a factor of 50 or more.

Further performance in your code would gained by restructuring the
computations so that harmonics up to a fixed order are computed at once,
as that is more efficient due to recursive definition.

ken-sakaie · 2014-08-14T17:35:37Z

Stefan

The code isn't copy and paste, so I think there shouldn't be copyright issues.

The algorithm really isn't much more than a recursion relation from Abramowitz and Stegun. A&S is now on the NIST website, which is US gov't funded (as was A&S). I expect it would be illegal to make the recursion relation proprietary.

Thank you,
Ken

On Aug 14, 2014, at 1:02 PM, "Stefan van der Walt" <notifications@github.com mailto:notifications@github.com> wrote:

Just note that numerical recipes code is not licensed under a free license
exactly, although I imagine it is intended to be used so.

On Thu, Aug 14, 2014 at 6:51 PM, ken-sakaie <notifications@github.com mailto:notifications@github.com>
wrote:

Eleftherios,

Couldn’t help but chime in on this one.

Numerical recipes has some really simple code for calculating these. See
the chapter on special functions, sub-chapter “Spherical Harmonics”

The code can be stripped down even more to speed things up.

I’ve attached some code I use (in C) to calculate the real spherical
harmonics.

Feel free to use it if it helps.

Sorry I don’t have time to do much coding.

Ken

From: Eleftherios Garyfallidis [mailto:notifications@github.com]
Sent: Wednesday, August 13, 2014 7:20 PM
To: nipy/dipy
Subject: Re: [dipy] A better way to create a response function for CSD
(#404)

@stefanvhttps://github.com/stefanv I remember that you had identified
in the past that real_sph_harm is too slow and you had some ideas on how to
better implement this. For @ChantalTaxhttps://github.com/ChantalTax the
performance of real_sph_harm is crucial as she is calling it many many
times. Would you mind suggesting a performance boost for this function.
Make it or show to Chantal how she could do it?

—
Reply to this email directly or view it on GitHub<
https://github.com/nipy/dipy/pull/404#issuecomment-52125000>.

Please consider the environment before printing this e-mail

Cleveland Clinic is ranked as one of the top hospitals in America by
U.S.News & World Report (2014).
Visit us online at http://www.clevelandclinic.org for a complete listing
of our services, staff and locations.

Confidentiality Note: This message is intended for use only by the
individual or entity to which it is addressed and may contain information
that is privileged, confidential, and exempt from disclosure under
applicable law. If the reader of this message is not the intended recipient
or the employee or agent responsible for delivering the message to the
intended recipient, you are hereby notified that any dissemination,
distribution or copying of this communication is strictly prohibited. If
you have received this communication in error, please contact the sender
immediately and destroy the material in its entirety, whether electronic or
hard copy.

Thank you.

/******************************************************************************
*

FILE: ylm_calc.c
*

USAGE: void ylm_calc(const int lmax, float z, float phi, float * ylm);
*

PURPOSE: Calculates real parts of spherical harmonic with antipodal

symmetry: Ylm(z=cos(theta),phi)
*

INPUT:
lmax: maximum degree, l. MUST BE EVEN.

z: cosine of angle with respect to z axis. ASSUMES -1<=z<=1

phi: azimuthal angle, IN RADIANS
*

OUTPUT:
ylm: array of spherical harmonics.

has (lmax+1)*(lmax/2+1) components.
*

Functions are in order:

N=1;

for L=0:LMAX

ylm(N) = YLM(L,0);

N=N+1;

for M=1:L

ylm(N) = -2*Im(YLM(L,M));

N=N+1;

ylm(N) = 2*Re(YLM(L,M));

N=N+1;

end

end
*
Copyright 2009 by Ken Sakaie. Please don't sue me.
$URL:
svn+ssh://pub1/home/sakaiek/repos/trunk/workspace/chardi/fod/ylm_calc.c $
$Author: sakaiek $
$Date: 2009-02-12 13:02:06 -0500 (Thu, 12 Feb 2009) $
$Rev: 1081 $
$Id: ylm_calc.c 1081 2009-02-12 18:02:06Z sakaiek $

*****************************************************************************/
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define M_4PI 12.56637061435917295405

void ylm_calc( int lmax, float z, float phi, float * ylm)
{
/* allocate workspace /
float * work = (float)malloc((lmax+1)*sizeof(float));
if(work == NULL) {
fprintf(stderr, "ylm_calc(work): memory error\n");
exit(1);
}

/* coefficient for legendre functions _/
int l;
for(l=0;l<=lmax;l++) work[l]=sqrt((2_l+1)/M_4PI);

/* variables /
int m;
float Pmm=1.0, mfac=1.0, lmfac=0.0; / for legendre /
float P_1=0.0, P_2=0.0, Pm1m=0.0, Plm=0.0; / for legendre recursion /
int j0=2,j=1; / indices for output /
int meven=1, leven=1; / flags: check m,l even /
float aylm=0.0, sm=0.0, cm=0.0; / function values */

/* calculate */
for(m=0;m<=lmax;m++) {

/* Ymm /
if ( m == 0 ) {
/ legendre /
ylm[0]=work[0];
/ set indices /
meven=0;
} else {
/ phase factors /
sm = -2.0 * sin( m * phi );
cm = 2.0 * cos( m * phi );
/ legendre _/
Pmm = -(2_m-1)_sqrt(1-z_z)Pmm;
mfac = mfac / sqrt( (m+m)(m+m-1) );
aylm = work[m] * mfac * Pmm;
/* index, assign _/
if ( meven==1 ) {
j = j0;
ylm[j] = aylm * sm;
ylm[j+1] = aylm * cm;
j0 += 2_m + 3;
j += 2*m + 1;
meven=0;
} else {
meven=1;
}
}

/* Ym+1,m /
if ( m < lmax ) {
/ legendre _/
Pm1m = z * (2_m+1) * Pmm;
lmfac = mfac / sqrt(m+m+1);
aylm = work[m+1] * lmfac * Pm1m;
/* index /
if ( meven==1 ) { / m odd, m+1 even /
j=j0;
ylm[j] = aylm * sm;
ylm[j+1] = aylm * cm;
j0 += 2;
j += 2(m+1)+1;
}
}

if ( m < (lmax-1) ) {
P_2 = Pmm;
P_1 = Pm1m;
leven = 1 - meven;
for(l=m+2;l<=lmax;l++) {
/* factorials /
lmfac *= sqrt(l-m) / sqrt(l+m);
/ legendre /
Plm = z( 2_l-1 )_P_1 - (l+m-1)P_2;
Plm = Plm / (l-m);
P_2 = P_1;
P_1 = Plm;
aylm = work[l] * lmfac * Plm;
/ index _/
if ( leven == 1 ) {
if ( m== 0 ) ylm[j] = aylm;
else {
ylm[j] = aylm * sm;
ylm[j+1] = aylm * cm;
}
j = j + 2_l + 1;
leven=0;
} else {
leven=1;
}
}
}
}

/* clean up */
free(work);

}
#undef M_4PI

—
Reply to this email directly or view it on GitHub
#404 (comment).

—

Reply to this email directly or view it on GitHubhttps://github.com//pull/404#issuecomment-52211303.

Please consider the environment before printing this e-mail

Cleveland Clinic is ranked as one of the top hospitals in America by U.S.News &
World Report (2014).
Visit us online at http://www.clevelandclinic.org for a complete listing of our
services, staff and locations.

Confidentiality Note: This message is intended for use only by the individual
or entity to which it is addressed and may contain information that is
privileged, confidential, and exempt from disclosure under applicable law. If
the reader of this message is not the intended recipient or the employee or
agent responsible for delivering the message to the intended recipient, you are
hereby notified that any dissemination, distribution or copying of this
communication is strictly prohibited. If you have received this communication
in error, please contact the sender immediately and destroy the material in its
entirety, whether electronic or hard copy.

Thank you.

pv · 2014-08-14T18:53:52Z

If you just want to calculate it for single m, n, use this (lpmv is pretty much the legendre recursion suggested above):

import numpy as np
from scipy.special import lpmv, gammaln

def sph_harm(m, n, theta, phi):
    x = np.cos(phi)
    val = lpmv(abs(m), n, x).astype(complex)
    val *= np.sqrt((2*n + 1) / 4.0 / np.pi)
    val *= np.exp(0.5*(gammaln(n-m+1)-gammaln(n+m+1)))
    val = val * np.exp(1j * m * theta)
    return val

The performance should be pretty close to C for large enough theta, phi arrays.

MrBago · 2014-08-21T18:57:01Z

dipy/reconst/csdeconv.py

+        peak in order to call it a single fiber population [1]
+    init_fa : float
+        FA of the initial 'fat' response function (tensor)
+    init_trace : fload


typo: fload -> float

MrBago · 2014-08-22T15:16:27Z

dipy/reconst/csdeconv.py

+    S0 = 1
+    evals = fa_trace_to_lambdas(init_fa, init_trace)
+    response = (evals, S0)
+    sphere = get_sphere('symmetric724')


To use a HemiSphere you can just replace this line with sphere = default_sphere and import default_sphere from dipy.reconst.peaks. default_sphere is HemiSphere version of symmetric724.

Please use this sphere from dipy.data import default_sphere.

MrBago · 2014-08-23T16:12:46Z

@ChantalTax if we add the faster implementation of the spherical harmonics as a new function, the rest of dipy will not benefit from the faster implementation (unless someone goes back and changes all the function calls). Can I suggest that if you'd like to work on a faster implementation, you split out this work into a separate pull request and simply replace sph_harm (or real_sph_harm if you have to). This pull request (response function) can be reviewed, tested and merged independent from how fast the spherical harmonic implementation is and I don't want the review of this very nice PR to digress into a conversion about how to implement spherical harmonic functions.

Also before you put a lot of time into this, I'm inclined to suggest that dipy continues to use the scipy implementation and that we make a note in the installation docs that recommends users upgrade to the upcoming scipy release for the best performance. I'm happy to discuss this point further but let's pick it up in issue #410.

arokem · 2014-10-06T16:42:08Z

Hey @ChantalTax - could you rebase this on master for the merge. I think this is otherwise ready to go, correct?

MrBago · 2014-10-07T02:42:40Z

I'd like to see the commits dealing with the shm basis rebased out of this before we merge it, but that shouldn't be much extra work because this needs to be rebased on master anyways.

arokem · 2014-10-07T16:07:02Z

Yeah - sounds OK to me. @ChantalTax, @Garyfallidis - please let us know if
you need any help rebasing this.

On Mon, Oct 6, 2014 at 7:42 PM, MrBago notifications@github.com wrote:

I'd like to see the commits dealing with the shm basis rebased out of this
before we merge it, but that shouldn't be much extra work because this
needs to be rebased on master anyways.

—
Reply to this email directly or view it on GitHub
#404 (comment).

Garyfallidis · 2014-10-07T19:10:12Z

@ChantalTax is currently out of office for the coming weeks. This PR will have to wait. I am sure she will need our help for rebasing this as she is new in the git business.

arokem · 2014-10-07T19:16:10Z

Should we just go ahead and rebase + merge this ourselves without her?

On Tue, Oct 7, 2014 at 12:10 PM, Eleftherios Garyfallidis <
notifications@github.com> wrote:

@ChantalTax https://github.com/ChantalTax is currently out of office
for the coming weeks. This PR will have to wait. I am sure she will need
our help for rebasing this as she is new in the git business.

—
Reply to this email directly or view it on GitHub
#404 (comment).

Garyfallidis · 2014-10-07T19:17:06Z

Better wait for Omar's PR to get in. This was rebased before his PR was reverted. It might resolve itself.

Garyfallidis · 2014-10-30T18:44:39Z

Hey, @ChantalTax are you back on this?

Garyfallidis · 2014-10-30T18:44:59Z

Let us know if you need any help with merging/rebasing etc.

MrBago · 2014-10-30T22:22:45Z

dipy/reconst/csdeconv.py

+
+            x, y, z = rot_gradients[where_dwi].T
+            r, theta, phi = cart2sphere(x, y, z)
+            # for the gradient sphere


I believe, but I might be wrong, that here you could do something like B_dwi = real_sph_harm(0, n[m == 0], theta[:, None], phi[:, None]). The coefficients associated with the other rows of the matrix are going to be set to zero anyways and the functions should be orthogonal, so the result should not change.

This would require changes in a few other places in the code, for example the CSD model might need accept the response function as N / 2 + 1 coefficients (N being the SH order) of a radially symmetric function. But that might be the right thing to do anyways.

Please let me know if I'm just wrong about this or if I'm not explaining it well.

kesshijordan · 2015-02-18T01:13:04Z

dipy/reconst/csdeconv.py

+        response[sh_mask] = 0
+
+        change = abs((response_p[~sh_mask] - response[~sh_mask])/response_p[~sh_mask])
+        if change.all() < convergence:


Hi, @ChantalTax

Is change.all() < convergence what you intended or did you mean all(change < convergence)? As it is now, the result of change.all() (either True or False) is compared to convergence (0.001 by default).

Cheers,
Kesshi

Hi @kesshijordan

Thanks for your comment, you are right, I think I mean the second option: to check if all changes were below the convergence threshold. I will replace this by all(change < convergence).

Best,
Chantal

arokem · 2015-03-11T16:44:08Z

Hey @ChantalTax - could you please rebase this on master? Looks like it won't merge automatically, so you might need to resolve some merge conflicts, before we can merge this. If you need any help with this, let me know - I'd be happy to walk you through this.

…ould be tested more because they don't give the same results as the currently used functions

… made for this

…esponse function axially symmetric

arokem · 2015-03-13T17:26:46Z

Hey @MrBago - this is now rebased on master. Could you please take another look?

arokem · 2015-03-13T18:43:18Z

dipy/reconst/csdeconv.py

@@ -47,7 +55,7 @@ def __init__(self, gtab, response, reg_sphere=None, sh_order=8, lambda_=1,
            ndarray and the second is the signal value for the response
            function without diffusion weighting.  This is to be able to
            generate a single fiber synthetic signal. The response function
-            will be used as deconvolution kernel ([1]_)
+            will be used as deconvolution git pull nipy-dipy masterkernel ([1]_)


Whoops. Looks like this somehow got inserted here, while we were rebasing this.

Garyfallidis · 2015-03-15T23:31:41Z

@ChantalTax, @arokem, @MrBago should we close this PR now?

Garyfallidis · 2015-03-15T23:32:17Z

Isn't that replaced by #588? Sorry, if I missed something.

MrBago · 2015-03-16T00:07:37Z

Thanks all! That's a wrap (We merged this in #588).

arokem reviewed Jul 31, 2014
View reviewed changes

MrBago reviewed Jul 31, 2014
View reviewed changes

MrBago reviewed Aug 21, 2014
View reviewed changes

MrBago mentioned this pull request Aug 21, 2014

Faster spherical harmonics implemenation #410

Closed

MrBago reviewed Aug 22, 2014
View reviewed changes

arokem mentioned this pull request Sep 6, 2014

Faster spherical harmonics #413

Merged

MrBago reviewed Oct 30, 2014
View reviewed changes

kesshijordan reviewed Feb 18, 2015
View reviewed changes

ChantalTax added 19 commits March 13, 2015 17:27

Starting working on csd calibration

22f5d0c

Added calculation of recursive response function

7199a4a

TEST: starting test on recursive calibration

72b916e

DOC: updating CSD example

3cc1890

RF: vectorized rotation of gradients

e598862

Fixed bug

e095d5d

Deleted visualization and fixed parallel processing

0787f84

TEST: testfile for recursive calibration

8c81afc

Removed not used lines

cc12b05

DOC: Extended manual with recursive calibration

9b91e63

added functions real_sph_harm2 and sph_harm2 which are faster, but sh…

555d171

…ould be tested more because they don't give the same results as the currently used functions

fixed bug on which Travis failed

e072a26

deleted faster implementation of sph_harm, a separate pull request is…

620ed14

… made for this

used hemisphere instead of sphere for peaks from model, made output r…

99e1a40

…esponse function axially symmetric

Updated with faster sph_harm version

b3df26a

Confirmed that new sph_harm is working much faster than before

30dddc0

fixed test, did not fail in approximately 100 trials

e807702

adapted convergence criterion

335e0e6

adding back code related to RF calibration

0c20dfd

ChantalTax force-pushed the csd_calibration branch from 2488ab2 to 0c20dfd Compare March 13, 2015 17:07

add back response scaling

bf7f292

arokem reviewed Mar 13, 2015
View reviewed changes

MrBago mentioned this pull request Mar 13, 2015

Nf csd calibration #588

Merged

MrBago closed this Mar 16, 2015

A better way to create a response function for CSD #404

A better way to create a response function for CSD #404

Conversation

ChantalTax commented Jul 31, 2014

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arokem commented Jul 31, 2014

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Garyfallidis commented Aug 13, 2014

stefanv commented Aug 13, 2014

Garyfallidis commented Aug 14, 2014

stefanv commented Aug 14, 2014

ken-sakaie commented Aug 14, 2014

stefanv commented Aug 14, 2014

pv commented Aug 14, 2014

ken-sakaie commented Aug 14, 2014

Reply to this email directly or view it on GitHubhttps://github.com//pull/404#issuecomment-52211303.

pv commented Aug 14, 2014

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

MrBago commented Aug 23, 2014

arokem commented Oct 6, 2014

MrBago commented Oct 7, 2014

arokem commented Oct 7, 2014

Garyfallidis commented Oct 7, 2014

arokem commented Oct 7, 2014

Garyfallidis commented Oct 7, 2014

Garyfallidis commented Oct 30, 2014

Garyfallidis commented Oct 30, 2014

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

arokem commented Mar 11, 2015

arokem commented Mar 13, 2015

Choose a reason for hiding this comment

Garyfallidis commented Mar 15, 2015

Garyfallidis commented Mar 15, 2015

MrBago commented Mar 16, 2015