bioem_algorithm.h 8.36 KB
Newer Older
Cossio, Pilar (pcos)'s avatar
Cossio, Pilar (pcos) committed
1
/* ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2
   < BioEM software for Bayesian inference of Electron Microscopy images>
Luka Stanisic's avatar
Luka Stanisic committed
3 4
   Copyright (C) 2017 Pilar Cossio, David Rohr, Fabio Baruffa, Markus Rampp,
        Luka Stanisic, Volker Lindenstruth and Gerhard Hummer.
5
   Max Planck Institute of Biophysics, Frankfurt, Germany.
Luka Stanisic's avatar
Luka Stanisic committed
6 7 8
   Frankfurt Institute for Advanced Studies, Goethe University Frankfurt,
   Germany.
   Max Planck Computing and Data Facility, Garching, Germany.
9

Luka Stanisic's avatar
Luka Stanisic committed
10
   Released under the GNU Public License, v3.
11
   See license statement for terms of distribution.
Cossio, Pilar (pcos)'s avatar
Cossio, Pilar (pcos) committed
12 13 14

   ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++*/

qon's avatar
qon committed
15 16
#ifndef BIOEM_ALGORITHM_H
#define BIOEM_ALGORITHM_H
17

Luka Stanisic's avatar
Luka Stanisic committed
18 19 20 21 22
__device__ static inline myprob_t
calc_logpro(const bioem_param_device &param, const myfloat_t amp,
            const myfloat_t pha, const myfloat_t env, const myfloat_t sum,
            const myfloat_t sumsquare, const myfloat_t crossproMapConv,
            const myfloat_t sumref, const myfloat_t sumsquareref)
23 24
{

25 26
  //*** MAIN ROUTINE TO CALCULATE THE LOGPRO FOR ALL KERNELS*************//
  // **** calculate the log posterior of Eq. of Pmw in SI of JSB paper ***//
27

28
  // Related to Reference calculated Projection
Luka Stanisic's avatar
Luka Stanisic committed
29
  const myprob_t ForLogProb = sumsquare * param.Ntotpi - sum * sum;
30

31
  // Products of different cross-correlations (first element in formula)
Luka Stanisic's avatar
Luka Stanisic committed
32 33 34 35 36
  const myprob_t firstele =
      param.Ntotpi *
          (sumsquareref * sumsquare - crossproMapConv * crossproMapConv) +
      2 * sumref * sum * crossproMapConv - sumsquareref * sum * sum -
      sumref * sumref * sumsquare;
37 38

  /// ******* Calculating log of Prob*********
Luka Stanisic's avatar
Luka Stanisic committed
39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
  // As in fortran code:
  // logpro=(3-Ntotpi)*0.5*log(firstele/pConvMap[iOrient].ForLogProbfromConv[iConv])+(Ntotpi*0.5-2)*log(Ntotpi-2)-0.5*log(pConvMap[iOrient].ForLogProbfromConv[iConv])+0.5*log(PI)+(1-Ntotpi*0.5)*(log(2*PI)+1);

  myprob_t logpro =
      (3 - param.Ntotpi) * 0.5 * log(firstele) +
      (param.Ntotpi * 0.5 - 2) * log((param.Ntotpi - 2) * ForLogProb);

  //*************Adding Gaussian Prior to envelope & Defocus
  // parameter******************

  if (not param.tousepsf)
  {
    logpro -= env * env / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf -
              (pha - param.Priordefcent) * (pha - param.Priordefcent) / 2. /
                  param.sigmaPriordefo / param.sigmaPriordefo -
              (amp - param.Priorampcent) * (amp - param.Priorampcent) / 2. /
                  param.sigmaPrioramp / param.sigmaPrioramp;
  }
  else
  {
    myprob_t envF, phaF;
    envF = 4. * M_PI * M_PI * env / (env * env + pha * pha);
    phaF = 4. * M_PI * M_PI * pha / (env * env + pha * pha);
    logpro -= envF * envF / 2. / param.sigmaPriorbctf / param.sigmaPriorbctf -
              (phaF - param.Priordefcent) * (phaF - param.Priordefcent) / 2. /
                  param.sigmaPriordefo / param.sigmaPriordefo -
              (amp - param.Priorampcent) * (amp - param.Priorampcent) / 2. /
                  param.sigmaPrioramp / param.sigmaPrioramp;
67
  }
68

Luka Stanisic's avatar
Luka Stanisic committed
69
  return (logpro);
70 71
}

Luka Stanisic's avatar
Luka Stanisic committed
72 73 74 75 76 77
__device__ static inline void
calProb(int iRefMap, int iOrient, int iConv, const myfloat_t amp,
        const myfloat_t pha, const myfloat_t env, const myfloat_t sumC,
        const myfloat_t sumsquareC, myfloat_t value, int disx, int disy,
        bioem_Probability &pProb, const bioem_param_device &param,
        const bioem_RefMap &RefMap)
78
{
79
  // IMPORTANT ROUTINE Summation of LogProb using FFTALGO
80 81 82
  // ********************************************************
  // *********** Calculates the BioEM probability ***********
  // ********************************************************
83

Luka Stanisic's avatar
Luka Stanisic committed
84 85 86
  myfloat_t logpro =
      calc_logpro(param, amp, pha, env, sumC, sumsquareC, value,
                  RefMap.sum_RefMap[iRefMap], RefMap.sumsquare_RefMap[iRefMap]);
87

Luka Stanisic's avatar
Luka Stanisic committed
88 89 90 91 92
#ifdef DEBUG_PROB
  printf("\t\t\tProb: iRefMap %d, iOrient %d, iConv %d, "
         "disx %d, disy %d, address -, value %f, logpro %f\n",
         iRefMap, iOrient, iConv, disx, disy, value, logpro);
#endif
93

Luka Stanisic's avatar
Luka Stanisic committed
94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122
  bioem_Probability_map &pProbMap = pProb.getProbMap(iRefMap);

  if (pProbMap.Constoadd < logpro)
  {
    pProbMap.Total *= exp(-logpro + pProbMap.Constoadd);
    pProbMap.Constoadd = logpro;

    // ********** Getting parameters that maximize the probability ***********
    pProbMap.max.max_prob_cent_x = -disx;
    pProbMap.max.max_prob_cent_y = -disy;
    pProbMap.max.max_prob_orient = iOrient;
    pProbMap.max.max_prob_conv = iConv;
    pProbMap.max.max_prob_norm =
        -(-sumC * RefMap.sum_RefMap[iRefMap] + param.Ntotpi * value) /
        (sumC * sumC - sumsquareC * param.Ntotpi);
    pProbMap.max.max_prob_mu =
        -(-sumC * value + sumsquareC * RefMap.sum_RefMap[iRefMap]) /
        (sumC * sumC - sumsquareC * param.Ntotpi);

#ifdef DEBUG_PROB
    printf("\tProbabilities change: iRefMap %d, iOrient %d, iConv %d, Total "
           "%f, Const %f, bestlogpro %f, sumExp -, bestId -\n",
           iRefMap, iOrient, iConv, pProbMap.Total, pProbMap.Constoadd, logpro);
    printf("\tParameters: iConv %d, myX -, myY -, disx %d, disy %d, probX "
           "%d, probY %d\n",
           iConv, disx, disy, pProbMap.max.max_prob_cent_x,
           pProbMap.max.max_prob_cent_y);
#endif
  }
123
  pProbMap.Total += exp(logpro - pProbMap.Constoadd);
Luka Stanisic's avatar
Luka Stanisic committed
124 125 126 127 128
#ifdef DEBUG_PROB
  printf("\t\tProbabilities after Sum: iRefMap %d, iOrient %d, iConv %d, "
         "Total %f, Const %f, bestlogpro %f, sumExp -, bestId -\n",
         iRefMap, iOrient, iConv, pProbMap.Total, pProbMap.Constoadd, logpro);
#endif
129 130

  if (param.writeAngles)
Luka Stanisic's avatar
Luka Stanisic committed
131 132
  {
    bioem_Probability_angle &pProbAngle = pProb.getProbAngle(iRefMap, iOrient);
133

Luka Stanisic's avatar
Luka Stanisic committed
134
    if (pProbAngle.ConstAngle < logpro)
135
    {
Luka Stanisic's avatar
Luka Stanisic committed
136 137
      pProbAngle.forAngles *= exp(-logpro + pProbAngle.ConstAngle);
      pProbAngle.ConstAngle = logpro;
138
    }
139

Luka Stanisic's avatar
Luka Stanisic committed
140 141
    pProbAngle.forAngles += exp(logpro - pProbAngle.ConstAngle);
  }
142 143
}

Luka Stanisic's avatar
Luka Stanisic committed
144 145 146 147 148 149
__device__ static inline void
doRefMapFFT(const int iRefMap, const int iOrient, const int iConv,
            const myfloat_t amp, const myfloat_t pha, const myfloat_t env,
            const myfloat_t sumC, const myfloat_t sumsquareC,
            const myfloat_t *lCC, bioem_Probability &pProb,
            const bioem_param_device &param, const bioem_RefMap &RefMap)
qon's avatar
qon committed
150
{
Luka Stanisic's avatar
Luka Stanisic committed
151 152 153 154 155 156 157 158 159 160
  //******************** Using FFT algorithm **************************
  //******************* Get cross-crollation of Ical to Iobs *******************
  //*********** Routine to get the Cross-Corellation from lCC for the interested
  // center displacement *************

  for (int cent_x = 0; cent_x <= param.maxDisplaceCenter;
       cent_x = cent_x + param.GridSpaceCenter)
  {
    for (int cent_y = 0; cent_y <= param.maxDisplaceCenter;
         cent_y = cent_y + param.GridSpaceCenter)
161
    {
Luka Stanisic's avatar
Luka Stanisic committed
162 163 164 165
      calProb(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC,
              (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] /
                  (myfloat_t)(param.NumberPixels * param.NumberPixels),
              cent_x, cent_y, pProb, param, RefMap);
166
    }
Luka Stanisic's avatar
Luka Stanisic committed
167 168
    for (int cent_y = param.NumberPixels - param.maxDisplaceCenter;
         cent_y < param.NumberPixels; cent_y = cent_y + param.GridSpaceCenter)
169
    {
Luka Stanisic's avatar
Luka Stanisic committed
170 171 172 173
      calProb(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC,
              (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] /
                  (myfloat_t)(param.NumberPixels * param.NumberPixels),
              cent_x, cent_y - param.NumberPixels, pProb, param, RefMap);
174
    }
Luka Stanisic's avatar
Luka Stanisic committed
175
  }
176

Luka Stanisic's avatar
Luka Stanisic committed
177 178 179 180 181
  for (int cent_x = param.NumberPixels - param.maxDisplaceCenter;
       cent_x < param.NumberPixels; cent_x = cent_x + param.GridSpaceCenter)
  {
    for (int cent_y = 0; cent_y <= param.maxDisplaceCenter;
         cent_y = cent_y + param.GridSpaceCenter)
182
    {
Luka Stanisic's avatar
Luka Stanisic committed
183 184 185 186
      calProb(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC,
              (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] /
                  (myfloat_t)(param.NumberPixels * param.NumberPixels),
              cent_x - param.NumberPixels, cent_y, pProb, param, RefMap);
187
    }
Luka Stanisic's avatar
Luka Stanisic committed
188 189
    for (int cent_y = param.NumberPixels - param.maxDisplaceCenter;
         cent_y < param.NumberPixels; cent_y = cent_y + param.GridSpaceCenter)
190
    {
Luka Stanisic's avatar
Luka Stanisic committed
191 192 193 194 195
      calProb(iRefMap, iOrient, iConv, amp, pha, env, sumC, sumsquareC,
              (myfloat_t) lCC[cent_x * param.NumberPixels + cent_y] /
                  (myfloat_t)(param.NumberPixels * param.NumberPixels),
              cent_x - param.NumberPixels, cent_y - param.NumberPixels, pProb,
              param, RefMap);
196
    }
Luka Stanisic's avatar
Luka Stanisic committed
197
  }
qon's avatar
qon committed
198 199 200
}

#endif