image_rotational_pca.cpp File Reference

#include "image_rotational_pca.h"
#include "core/metadata_extension.h"
#include "core/transformations.h"
#include "core/xmipp_funcs.h"

Include dependency graph for image_rotational_pca.cpp:

Go to the source code of this file.

Functions
void	threadApplyT (ThreadArgument &thArg)
void	threadApplyTt (ThreadArgument &thArg)

Function Documentation

threadApplyT()

void threadApplyT

(

ThreadArgument &

thArg

)

Definition at line 227 of file image_rotational_pca.cpp.

{
  auto *self=(ProgImageRotationalPCA *) thArg.workClass;
  //MpiNode *node=self->node;
  int rank = self->rank;
  std::vector<size_t> &objId=self->objId;
  MetaDataVec &MD=self->MD[thArg.thread_id];

  Image<double> &I=self->I[thArg.thread_id];
  MultidimArray<double> &Iaux=self->Iaux[thArg.thread_id];
  Matrix2D<double> &Wnode=self->Wnode[thArg.thread_id];
  Matrix2D<double> &Hblock=self->Hblock[thArg.thread_id];
  Matrix2D<double> &A=self->A[thArg.thread_id];
  Matrix2D<double> &H=self->H;
  MultidimArray< unsigned char > &mask=self->mask;
  Wnode.initZeros(self->Npixels,MAT_XSIZE(H));

  const int unroll=8;
  const int jmax=(MAT_XSIZE(Wnode)/unroll)*unroll;

  size_t first, last;
  if (IS_MASTER && thArg.thread_id==0)
  {
    std::cout << "Applying T ...\n";
    init_progress_bar(objId.size());
  }
  while (self->taskDistributor->getTasks(first, last))
  {
    for (size_t idx=first; idx<=last; ++idx)
    {
      // Read image
      I.readApplyGeo(MD,objId[idx]);
      MultidimArray<double> &mI=I();

      // Locate the corresponding index in Matrix H
      // and copy a block in memory to speed up calculations
      size_t Hidx=idx*2*self->Nangles*self->Nshifts;
      memcpy(&MAT_ELEM(Hblock,0,0),&MAT_ELEM(H,Hidx,0),MAT_XSIZE(Hblock)*MAT_YSIZE(Hblock)*sizeof(double));

      // For each rotation, shift and mirror
      int block_idx=0;
      for (int mirror=0; mirror<2; ++mirror)
      {
        if (mirror)
        {
          mI.selfReverseX();
          mI.setXmippOrigin();
        }
        for (double psi=0; psi<360; psi+=self->psi_step)
        {
          rotation2DMatrix(psi,A,true);
          for (double y=-self->max_shift_change; y<=self->max_shift_change; y+=self->shift_step)
          {
            MAT_ELEM(A,1,2)=y;
            for (double x=-self->max_shift_change; x<=self->max_shift_change; x+=self->shift_step, ++block_idx)
            {
              MAT_ELEM(A,0,2)=x;

              // Rotate and shift image
              applyGeometry(xmipp_transformation::LINEAR,Iaux,mI,A,xmipp_transformation::IS_INV,true);

              // Update Wnode
              int i=0;
              FOR_ALL_DIRECT_ELEMENTS_IN_MULTIDIMARRAY(Iaux)
              {
                if (DIRECT_MULTIDIM_ELEM(mask,n)==0)
                continue;
                double pixval=DIRECT_MULTIDIM_ELEM(Iaux,n);
                double *ptrWnode=&MAT_ELEM(Wnode,i,0);
                double *ptrHblock=&MAT_ELEM(Hblock,block_idx,0);
                for (int j=0; j<jmax; j+=unroll, ptrHblock+=unroll, ptrWnode+=unroll)
                {
                  (*ptrWnode) +=pixval*(*ptrHblock);
                  (*(ptrWnode+1)) +=pixval*(*(ptrHblock+1));
                  (*(ptrWnode+2)) +=pixval*(*(ptrHblock+2));
                  (*(ptrWnode+3)) +=pixval*(*(ptrHblock+3));
                  (*(ptrWnode+4)) +=pixval*(*(ptrHblock+4));
                  (*(ptrWnode+5)) +=pixval*(*(ptrHblock+5));
                  (*(ptrWnode+6)) +=pixval*(*(ptrHblock+6));
                  (*(ptrWnode+7)) +=pixval*(*(ptrHblock+7));
                }
                for (size_t j=jmax; j<MAT_XSIZE(Wnode); ++j, ptrHblock+=1, ptrWnode+=1)
                (*ptrWnode) +=pixval*(*ptrHblock );
                ++i;
              }
            }
          }
        }
      }
    }
    if (IS_MASTER && thArg.thread_id==0)
    progress_bar(last);
  }
}

threadApplyTt()

void threadApplyTt

(

ThreadArgument &

thArg

)

Definition at line 342 of file image_rotational_pca.cpp.

{
  auto *self=(ProgImageRotationalPCA *) thArg.workClass;
  //MpiNode *node=self->node;
  int rank = self->rank;
  std::vector<size_t> &objId=self->objId;
  MetaDataVec &MD=self->MD[thArg.thread_id];

  Image<double> &I=self->I[thArg.thread_id];
  MultidimArray<double> &Iaux=self->Iaux[thArg.thread_id];
  Matrix2D<double> &A=self->A[thArg.thread_id];
  Matrix2D<double> &Hblock=self->Hblock[thArg.thread_id];
  Matrix2D<double> &H=self->H;
  Matrix2D<double> &Wtranspose=self->Wtranspose;
  MultidimArray< unsigned char > &mask=self->mask;

  const size_t unroll=8;
  const size_t nmax=(MAT_XSIZE(Wtranspose)/unroll)*unroll;

  size_t first, last;
  if (IS_MASTER && thArg.thread_id==0)
  {
    std::cout << "Applying Tt ...\n";
    init_progress_bar(objId.size());
  }
  while (self->taskDistributor->getTasks(first, last))
  {
    for (size_t idx=first; idx<=last; ++idx)
    {
      // Read image
      I.readApplyGeo(MD,objId[idx]);
      MultidimArray<double> &mI=I();

      // For each rotation and shift
      int block_idx=0;
      for (int mirror=0; mirror<2; ++mirror)
      {
        if (mirror)
        {
          mI.selfReverseX();
          mI.setXmippOrigin();
        }
        for (double psi=0; psi<360; psi+=self->psi_step)
        {
          rotation2DMatrix(psi,A,true);
          for (double y=-self->max_shift_change; y<=self->max_shift_change; y+=self->shift_step)
          {
            MAT_ELEM(A,1,2)=y;
            for (double x=-self->max_shift_change; x<=self->max_shift_change; x+=self->shift_step, ++block_idx)
            {
              MAT_ELEM(A,0,2)=x;

              // Rotate and shift image
              applyGeometry(xmipp_transformation::LINEAR,Iaux,mI,A,xmipp_transformation::IS_INV,true);

              // Update Hblock
              for (size_t j=0; j<MAT_XSIZE(Hblock); j++)
              {
                double dotproduct=0;
                const double *ptrIaux=MULTIDIM_ARRAY(Iaux);
                unsigned char *ptrMask=&DIRECT_MULTIDIM_ELEM(mask,0);
                const double *ptrWtranspose=&MAT_ELEM(Wtranspose,j,0);
                for (size_t n=0; n<nmax; n+=unroll, ptrIaux+=unroll, ptrMask+=unroll)
                {
                  if (*ptrMask)
                  dotproduct+=(*ptrIaux)*(*ptrWtranspose++);
                  if (*(ptrMask+1))
                  dotproduct+=(*(ptrIaux+1))*(*ptrWtranspose++);
                  if (*(ptrMask+2))
                  dotproduct+=(*(ptrIaux+2))*(*ptrWtranspose++);
                  if (*(ptrMask+3))
                  dotproduct+=(*(ptrIaux+3))*(*ptrWtranspose++);
                  if (*(ptrMask+4))
                  dotproduct+=(*(ptrIaux+4))*(*ptrWtranspose++);
                  if (*(ptrMask+5))
                  dotproduct+=(*(ptrIaux+5))*(*ptrWtranspose++);
                  if (*(ptrMask+6))
                  dotproduct+=(*(ptrIaux+6))*(*ptrWtranspose++);
                  if (*(ptrMask+7))
                  dotproduct+=(*(ptrIaux+7))*(*ptrWtranspose++);
                }
                for (size_t n=nmax; n<MAT_XSIZE(Wtranspose); ++n, ++ptrMask, ++ptrIaux)
                if (*ptrMask)
                dotproduct+=(*ptrIaux)*(*ptrWtranspose++);
                MAT_ELEM(Hblock,block_idx,j)=dotproduct;
              }
            }
          }
        }
      }

      // Locate the corresponding index in Matrix H
      // and copy block to disk
      size_t Hidx=idx*2*self->Nangles*self->Nshifts;
      // COSS std::cout << "idx=" << idx << " Hidx=" << Hidx << std::endl;
      // COSS std::cout << "H shape " << MAT_YSIZE(H) << " x " << MAT_XSIZE(H) << std::endl;
      self->writeToHBuffer(thArg.thread_id,&MAT_ELEM(H,Hidx,0));
    }
    if (IS_MASTER && thArg.thread_id==0)
      progress_bar(last);
  }
}