Xmipp  v3.23.11-Nereus
Macros | Functions | Variables
blobs.cpp File Reference
#include <pthread.h>
#include "blobs.h"
#include "core/numerical_recipes.h"
#include "xmipp_image_over.h"
Include dependency graph for blobs.cpp:

Go to the source code of this file.

Macros

#define DEFORM_BLOB_WHEN_IN_CRYSTAL
 
#define x0   STARTINGX(*vol_voxels)
 
#define xF   FINISHINGX(*vol_voxels)
 
#define y0   STARTINGY(*vol_voxels)
 
#define yF   FINISHINGY(*vol_voxels)
 
#define z0   STARTINGZ(*vol_voxels)
 
#define zF   FINISHINGZ(*vol_voxels)
 
#define FORWARD   true
 
#define BACKWARD   false
 

Functions

double kaiser_value (double r, double a, double alpha, int m)
 
double kaiser_proj (double s, double a, double alpha, int m)
 
double kaiser_Fourier_value (double w, double a, double alpha, int m)
 
double sum_blob_SimpleGrid (const struct blobtype &blob, const SimpleGrid &grid, const Matrix2D< double > *D)
 
double basvolume (double a, double alpha, int m, int n)
 
double i_n (int n, double x)
 
double i_nph (int n, double x)
 
double in_zeroarg (int n)
 
double inph_zeroarg (int n)
 
double sum_blob_Grid (const struct blobtype &blob, const Grid &grid, const Matrix2D< double > *D)
 
double blob_freq_zero (struct blobtype b)
 
double blob_att (double w, struct blobtype b)
 
double blob_ops (double w, struct blobtype b)
 
double optimal_CC_grid_relative_size (struct blobtype b)
 
double optimal_BCC_grid_relative_size (struct blobtype b)
 
double optimal_FCC_grid_relative_size (struct blobtype b)
 
blobtype best_blob (double alpha_0, double alpha_F, double inc_alpha, double a_0, double a_F, double inc_a, double w, double *target_att, int target_length)
 
void footprint_blob (ImageOver &blobprint, const struct blobtype &blob, int istep, int normalise)
 
void * blobs2voxels_SimpleGrid (void *data)
 
void voxel_volume_shape (const GridVolume &vol_blobs, const struct blobtype &blob, const Matrix2D< double > *D, Matrix1D< int > &corner1, Matrix1D< int > &size)
 
void blobs2voxels (const GridVolume &vol_blobs, const struct blobtype &blob, MultidimArray< double > *vol_voxels, const Matrix2D< double > *D, int threads, int Zdim, int Ydim, int Xdim)
 
void blobs2space_coefficients (const GridVolume &vol_blobs, const struct blobtype &blob, MultidimArray< double > *vol_coefs)
 
void ART_voxels2blobs_single_step (GridVolume &vol_in, GridVolume *vol_out, const struct blobtype &blob, const Matrix2D< double > *D, double lambda, MultidimArray< double > *theo_vol, const MultidimArray< double > *read_vol, MultidimArray< double > *corr_vol, const MultidimArray< double > *mask_vol, double &mean_error, double &max_error, int eq_mode, int threads)
 
void voxels2blobs (const MultidimArray< double > *vol_voxels, const struct blobtype &blob, GridVolume &vol_blobs, int grid_type, double grid_relative_size, double lambda, const MultidimArray< double > *vol_mask, const Matrix2D< double > *D, double final_error_change, int tell, double R, int threads)
 

Variables

pthread_mutex_t blobs_conv_mutex = PTHREAD_MUTEX_INITIALIZER
 
int * slices_status
 
int slices_processed
 

Macro Definition Documentation

◆ BACKWARD

#define BACKWARD   false

Definition at line 1092 of file blobs.cpp.

◆ DEFORM_BLOB_WHEN_IN_CRYSTAL

#define DEFORM_BLOB_WHEN_IN_CRYSTAL

Definition at line 451 of file blobs.cpp.

◆ FORWARD

#define FORWARD   true

Definition at line 1091 of file blobs.cpp.

◆ x0

#define x0   STARTINGX(*vol_voxels)

◆ xF

#define xF   FINISHINGX(*vol_voxels)

◆ y0

#define y0   STARTINGY(*vol_voxels)

◆ yF

#define yF   FINISHINGY(*vol_voxels)

◆ z0

#define z0   STARTINGZ(*vol_voxels)

◆ zF

#define zF   FINISHINGZ(*vol_voxels)

Function Documentation

◆ blobs2voxels_SimpleGrid()

void* blobs2voxels_SimpleGrid ( void *  data)

Definition at line 452 of file blobs.cpp.

453 {
454  auto * thread_data = (ThreadBlobsToVoxels *) data;
455 
456  const MultidimArray<double> *vol_blobs = thread_data->vol_blobs;
457  const SimpleGrid *grid = thread_data->grid;
458  const struct blobtype *blob = thread_data->blob;
459  MultidimArray<double> *vol_voxels = thread_data->vol_voxels;
460  const Matrix2D<double> *D = thread_data->D;
461  int istep = thread_data->istep;
462  MultidimArray<double> *vol_corr = thread_data->vol_corr;
463  const MultidimArray<double> *vol_mask = thread_data->vol_mask;
464  ;
465  bool FORW = thread_data->FORW;
466  int eq_mode = thread_data->eq_mode;
467 
468  int min_separation = thread_data->min_separation;
469 
470  auto z_planes = (int)(ZZ(grid->highest) - ZZ(grid->lowest) + 1);
471 
472  Matrix2D<double> Dinv; // Inverse of D
473  Matrix1D<double> act_coord(3); // Coord: Actual position inside
474  // the voxel volume without deforming
475  Matrix1D<double> real_position(3); // Coord: actual position after
476  // applying the V transformation
477  Matrix1D<double> beginZ(3); // Coord: Voxel coordinates of the
478  // blob at the 3D point
479  // (z0,YY(lowest),XX(lowest))
480  Matrix1D<double> beginY(3); // Coord: Voxel coordinates of the
481  // blob at the 3D point
482  // (z0,y0,XX(lowest))
483  Matrix1D<double> corner2(3);
484  Matrix1D<double> corner1(3); // Coord: Corners of the blob in the voxel volume
485  Matrix1D<double> gcurrent(3); // Position in g of current point
486  MultidimArray<double> blob_table; // Something like a blobprint but with the values of the
487  // blob in space
488  double d; // Distance between the center
489  // of the blob and a voxel position
490  int id; // index inside the blob value
491  // table for the blob value at a distance d
492  int i;
493  int j;
494  int k; // Index within the blob volume
495  int process; // True if this blob has to be
496  // processed
497  double vol_correction=0; // Correction to apply to the
498  // volume when "projecting" back
500 
501  // Some aliases
502 #define x0 STARTINGX(*vol_voxels)
503 #define xF FINISHINGX(*vol_voxels)
504 #define y0 STARTINGY(*vol_voxels)
505 #define yF FINISHINGY(*vol_voxels)
506 #define z0 STARTINGZ(*vol_voxels)
507 #define zF FINISHINGZ(*vol_voxels)
508 
509 #ifdef DEBUG
510 
511  bool condition = !FORW;
512  if (condition)
513  {
514  (*vol_voxels)().printShape();
515  std::cout << std::endl;
516  std::cout << "x0= " << x0 << " xF= " << xF << std::endl;
517  std::cout << "y0= " << y0 << " yF= " << yF << std::endl;
518  std::cout << "z0= " << z0 << " zF= " << zF << std::endl;
519  std::cout << grid;
520  }
521 #endif
522 
523  // Invert deformation matrix ............................................
524  if (D != nullptr)
525  Dinv = D->inv();
526 
527  // Compute a blob value table ...........................................
528  blob_table.resize((int)(blob->radius*istep + 1));
529  for (size_t i = 0; i < blob_table.xdim; i++)
530  {
531  A1D_ELEM(blob_table, i) = kaiser_value((double)i/istep, blob->radius, blob->alpha, blob->order);
532 
533 #ifdef DEBUG_MORE
534 
535  if (condition)
536  std::cout << "Blob (" << i << ") r=" << (double)i / istep <<
537  " val= " << A1D_ELEM(blob_table, i) << std::endl;
538 #endif
539 
540  }
541 
542  int assigned_slice;
543 
544  do
545  {
546  assigned_slice = -1;
547  do
548  {
549  pthread_mutex_lock(&blobs_conv_mutex);
550  if( slices_processed == z_planes )
551  {
552  pthread_mutex_unlock(&blobs_conv_mutex);
553  return (void*)nullptr;
554  }
555 
556  for(int w = 0 ; w < z_planes ; w++ )
557  {
558  if( slices_status[w]==0 )
559  {
560  slices_status[w] = -1;
561  assigned_slice = w;
563 
564  for( int in = (w-min_separation+1) ; in <= (w+min_separation-1 ) ; in ++ )
565  {
566  if( in != w )
567  {
568  if( ( in >= 0 ) && ( in < z_planes ))
569  {
570  if( slices_status[in] != -1 )
571  slices_status[in]++;
572  }
573  }
574  }
575  break;
576  }
577  }
578 
579  pthread_mutex_unlock(&blobs_conv_mutex);
580  }
581  while( assigned_slice == -1);
582 
583  // Convert the whole grid ...............................................
584  // Corner of the plane defined by Z. These coordinates are in the
585  // universal coord. system
586  Matrix1D<double> aux( grid->lowest );
587  k = (int)(assigned_slice + ZZ( grid->lowest ));
588  ZZ(aux) = k;
589  grid->grid2universe(aux, beginZ);
590 
591  Matrix1D<double> grid_index(3);
592 
593  // Corner of the row defined by Y
594  beginY = beginZ;
595  for (i = (int) YY(grid->lowest); i <= (int) YY(grid->highest); i++)
596  {
597  // First point in the row
598  act_coord = beginY;
599  for (j = (int) XX(grid->lowest); j <= (int) XX(grid->highest); j++)
600  {
601  VECTOR_R3(grid_index, j, i, k);
602 #ifdef DEBUG
603 
604  if (condition)
605  {
606  printf("Dealing blob at (%d,%d,%d) = %f\n", j, i, k, A3D_ELEM(*vol_blobs, k, i, j));
607  std::cout << "Center of the blob "
608  << act_coord.transpose() << std::endl;
609  }
610 #endif
611 
612  // Place act_coord in its right place
613  if (D != nullptr)
614  {
615  M3x3_BY_V3x1(real_position, *D, act_coord);
616 #ifdef DEBUG
617 
618  if (condition)
619  std::cout << "Center of the blob moved to "
620  //ROB, the "moved" coordinates are in
621  // real_position not in act_coord
622  << act_coord.transpose() << std::endl;
623  << real_position.transpose() << std::endl;
624 #endif
625  // ROB This is OK if blob.radius is in Cartesian space as I
626  // think is the case
627  }
628  else
629  real_position = act_coord;
630 
631  // These two corners are also real valued
632  process = true;
633  //ROB
634  //This is OK if blob.radius is in Cartesian space as I think is the case
635  V3_PLUS_CT(corner1, real_position, -blob->radius);
636  V3_PLUS_CT(corner2, real_position, blob->radius);
637 //#ifdef DEFORM_BLOB_WHEN_IN_CRYSTAL
638  //ROB
639  //we do not need this, it is already in Cartesian space
640  //if (D!=NULL)
641  // box_enclosing(corner1,corner2, *D, corner1, corner2);
642 //#endif
643 
644  if (XX(corner1) >= xF)
645  process = false;
646  if (YY(corner1) >= yF)
647  process = false;
648  if (ZZ(corner1) >= zF)
649  process = false;
650  if (XX(corner2) <= x0)
651  process = false;
652  if (YY(corner2) <= y0)
653  process = false;
654  if (ZZ(corner2) <= z0)
655  process = false;
656 #ifdef DEBUG
657 
658  if (!process && condition)
659  std::cout << " It is outside output volume\n";
660 #endif
661 
662  if (!grid->is_interesting(real_position))
663  {
664 #ifdef DEBUG
665  if (process && condition)
666  std::cout << " It is not interesting\n";
667 #endif
668 
669  process = false;
670  }
671 
672 #ifdef DEBUG
673  if (condition)
674  {
675  std::cout << "Corner 1 for this point " << corner1.transpose() << std::endl;
676  std::cout << "Corner 2 for this point " << corner2.transpose() << std::endl;
677  }
678 #endif
679 
680  if (process)
681  {
682  // Clip the corners to the volume borders
683  XX(corner1) = ROUND(CLIP(XX(corner1), x0, xF));
684  YY(corner1) = ROUND(CLIP(YY(corner1), y0, yF));
685  ZZ(corner1) = ROUND(CLIP(ZZ(corner1), z0, zF));
686  XX(corner2) = ROUND(CLIP(XX(corner2), x0, xF));
687  YY(corner2) = ROUND(CLIP(YY(corner2), y0, yF));
688  ZZ(corner2) = ROUND(CLIP(ZZ(corner2), z0, zF));
689 #ifdef DEBUG
690 
691  if (condition)
692  {
693  std::cout << "Clipped and rounded Corner 1 " << corner1.transpose() << std::endl;
694  std::cout << "Clipped and rounded Corner 2 " << corner2.transpose() << std::endl;
695  }
696 #endif
697 
698  if (!FORW)
699  switch (eq_mode)
700  {
701  case VARTK:
702  vol_correction = 0;
703  break;
704  case VMAXARTK:
705  vol_correction = -1e38;
706  break;
707  }
708 
709  // Effectively convert
710  long N_eq;
711  N_eq = 0;
712  for (auto intz = (int)ZZ(corner1); intz <=(int)ZZ(corner2); intz++)
713  for (auto inty = (int)YY(corner1); inty <= (int)YY(corner2); inty++)
714  for (auto intx = (int)XX(corner1); intx <= (int)XX(corner2); intx++)
715  {
716  if (vol_mask != nullptr && A3D_ELEM(*vol_mask, intz, inty, intx)!=0.0)
717  continue;
718 
719  // Compute distance to the center of the blob
720  VECTOR_R3(gcurrent, (double)intx, (double)inty, (double)intz);
721 //#ifdef DEFORM_BLOB_WHEN_IN_CRYSTAL
722  // ROB
723  //if (D!=NULL)
724  // M3x3_BY_V3x1(gcurrent,Dinv,gcurrent);
725 //#endif
726 
727  V3_MINUS_V3(gcurrent, real_position, gcurrent);
728  d = sqrt(XX(gcurrent) * XX(gcurrent) +
729  YY(gcurrent) * YY(gcurrent) +
730  ZZ(gcurrent) * ZZ(gcurrent));
731  if (d > blob->radius)
732  continue;
733  id = (int)(d * istep);
734 #ifdef DEBUG_MORE
735 
736  if (condition)
737  {
738  std::cout << "At (" << intx << ","
739  << inty << "," << intz << ") distance=" << d;
740  std::cout.flush();
741  }
742 #endif
743 
744  // Add at that position the corresponding blob value
745 
746  if (FORW)
747  {
748  A3D_ELEM(*vol_voxels, intz, inty, intx) +=
749  A3D_ELEM(*vol_blobs, k, i, j) *
750  A1D_ELEM(blob_table, id);
751 #ifdef DEBUG_MORE
752 
753  if (condition)
754  {
755  std::cout << " adding " << A3D_ELEM(*vol_blobs, k, i, j)
756  << " * " << A1D_ELEM(blob_table, id) << " = "
757  << A3D_ELEM(*vol_blobs, k, i, j)*
758  A1D_ELEM(blob_table, id) << std::endl;
759  std::cout.flush();
760  }
761 #endif
762  if (vol_corr != nullptr)
763  A3D_ELEM(*vol_corr, intz, inty, intx) +=
764  A1D_ELEM(blob_table, id) * A1D_ELEM(blob_table, id);
765  }
766  else
767  {
768  double contrib = A3D_ELEM(*vol_corr, intz, inty, intx) *
769  A1D_ELEM(blob_table, id);
770  switch (eq_mode)
771  {
772  case VARTK:
773  vol_correction += contrib;
774  N_eq++;
775  break;
776  case VMAXARTK:
777  if (contrib > vol_correction)
778  vol_correction = contrib;
779  break;
780 
781  }
782 #ifdef DEBUG_MORE
783  if (condition)
784  {
785  std::cout << " adding " << A3D_ELEM(*vol_corr, intz, inty, intx)
786  << " * " << A1D_ELEM(blob_table, id) << " = "
787  << contrib << std::endl;
788  std::cout.flush();
789  }
790 #endif
791 
792  }
793  }
794  if (N_eq == 0)
795  N_eq = 1;
796  if (!FORW)
797  {
798  A3D_ELEM(*vol_blobs, k, i, j) += vol_correction / N_eq;
799 #ifdef DEBUG_MORE
800 
801  std::cout << " correction= " << vol_correction << std::endl
802  << " Number of eqs= " << N_eq << std::endl
803  << " Blob after correction= "
804  << A3D_ELEM(*vol_blobs, k, i, j) << std::endl;
805 #endif
806 
807  }
808  }
809 
810  // Prepare for next iteration
811  XX(act_coord) = XX(act_coord) + grid->relative_size * (grid->basis)( 0, 0);
812  YY(act_coord) = YY(act_coord) + grid->relative_size * (grid->basis)( 1, 0);
813  ZZ(act_coord) = ZZ(act_coord) + grid->relative_size * (grid->basis)( 2, 0);
814  }
815  XX(beginY) = XX(beginY) + grid->relative_size * (grid->basis)( 0, 1);
816  YY(beginY) = YY(beginY) + grid->relative_size * (grid->basis)( 1, 1);
817  ZZ(beginY) = ZZ(beginY) + grid->relative_size * (grid->basis)( 2, 1);
818  }
819 
820  pthread_mutex_lock(&blobs_conv_mutex);
821 
822  for( int in = (assigned_slice-min_separation+1) ; in <= (assigned_slice+min_separation-1); in ++ )
823  {
824  if( in != assigned_slice )
825  {
826  if( ( in >= 0 ) && ( in < z_planes))
827  {
828  if( slices_status[in] != -1 )
829  slices_status[in]--;
830  }
831  }
832  }
833 
834  pthread_mutex_unlock(&blobs_conv_mutex);
835  }
836  while(1);
837 
838  return (void*)nullptr;
839 }
double kaiser_value(double r, double a, double alpha, int m)
Definition: blobs.cpp:37
double alpha
Smoothness parameter.
Definition: blobs.h:121
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
Matrix1D< double > highest
Definition: grids.h:161
#define V3_PLUS_CT(a, b, c)
Definition: matrix1d.h:220
void grid2universe(const Matrix1D< double > &gv, Matrix1D< double > &uv) const
Definition: grids.h:318
void sqrt(Image< double > &op)
constexpr int VARTK
Definition: blobs.h:378
#define V3_MINUS_V3(a, b, c)
Definition: matrix1d.h:202
#define A1D_ELEM(v, i)
pthread_mutex_t blobs_conv_mutex
Definition: blobs.cpp:31
void inv(Matrix2D< T > &result) const
Definition: matrix2d.cpp:663
doublereal * w
#define zF
#define z0
bool is_interesting(const Matrix1D< double > &uv) const
Definition: grids.h:361
Matrix1D< double > lowest
Definition: grids.h:158
Matrix2D< double > basis
Definition: grids.h:148
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
doublereal * d
int * slices_status
Definition: blobs.cpp:33
#define CLIP(x, x0, xF)
Definition: xmipp_macros.h:260
#define A3D_ELEM(V, k, i, j)
#define y0
#define x0
#define XX(v)
Definition: matrix1d.h:85
int in
#define ROUND(x)
Definition: xmipp_macros.h:210
#define xF
double relative_size
Measuring unit in the grid coordinate system.
Definition: grids.h:163
#define j
#define M3x3_BY_V3x1(a, M, b)
Definition: matrix2d.h:170
#define YY(v)
Definition: matrix1d.h:93
#define VECTOR_R3(v, x, y, z)
Definition: matrix1d.h:124
int order
Derivation order and Bessel function order.
Definition: blobs.h:118
int slices_processed
Definition: blobs.cpp:34
#define yF
double radius
Spatial radius in Universal System units.
Definition: blobs.h:115
#define ZZ(v)
Definition: matrix1d.h:101
constexpr int VMAXARTK
Definition: blobs.h:379
#define SPEED_UP_temps012
Definition: xmipp_macros.h:403

◆ sum_blob_SimpleGrid()

double sum_blob_SimpleGrid ( const struct blobtype blob,
const SimpleGrid grid,
const Matrix2D< double > *  D 
)

Definition at line 175 of file blobs.cpp.

177 {
179  Matrix1D<double> gr(3);
180  Matrix1D<double> ur(3);
181  Matrix1D<double> corner1(3);
182  Matrix1D<double> corner2(3);
183  double actual_radius;
184  int i;
185  int j;
186  int k;
187  double sum = 0.0;
188 
189  // Compute the limits of the blob in the grid coordinate system
190  grid.universe2grid(vectorR3(-blob.radius, -blob.radius, -blob.radius), corner1);
191  grid.universe2grid(vectorR3(blob.radius, blob.radius, blob.radius), corner2);
192  if (D != nullptr)
193  box_enclosing(corner1, corner2, *D, corner1, corner2);
194 
195  // Compute the sum in the points inside the grid
196  // The integer part of the vectors is taken for not picking points
197  // just in the border of the blob, which we know they are 0.
198  for (i = (int)XX(corner1); i <= (int)XX(corner2); i++)
199  for (j = (int)YY(corner1); j <= (int)YY(corner2); j++)
200  for (k = (int)ZZ(corner1); k <= (int)ZZ(corner2); k++)
201  {
202  VECTOR_R3(gr, i, j, k);
203  grid.grid2universe(gr, ur);
204  if (D != nullptr)
205  M3x3_BY_V3x1(ur, *D, ur);
206  actual_radius = ur.module();
207  if (actual_radius < blob.radius)
208  sum += kaiser_value(actual_radius,
209  blob.radius, blob.alpha, blob.order);
210  }
211  return sum;
212 }
double kaiser_value(double r, double a, double alpha, int m)
Definition: blobs.cpp:37
void universe2grid(const Matrix1D< double > &uv, Matrix1D< double > &gv) const
Definition: grids.h:349
double alpha
Smoothness parameter.
Definition: blobs.h:121
void grid2universe(const Matrix1D< double > &gv, Matrix1D< double > &uv) const
Definition: grids.h:318
Matrix1D< double > vectorR3(double x, double y, double z)
Definition: matrix1d.cpp:892
#define i
ql0001_ & k(htemp+1),(cvec+1),(atemp+1),(bj+1),(bl+1),(bu+1),(x+1),(clamda+1), &iout, infoqp, &zero,(w+1), &lenw,(iw+1), &leniw, &glob_grd.epsmac
void box_enclosing(const Matrix1D< double > &v0, const Matrix1D< double > &vF, const Matrix2D< double > &V, Matrix1D< double > &corner1, Matrix1D< double > &corner2)
Definition: geometry.cpp:366
#define XX(v)
Definition: matrix1d.h:85
#define j
#define M3x3_BY_V3x1(a, M, b)
Definition: matrix2d.h:170
#define YY(v)
Definition: matrix1d.h:93
#define VECTOR_R3(v, x, y, z)
Definition: matrix1d.h:124
int order
Derivation order and Bessel function order.
Definition: blobs.h:118
double radius
Spatial radius in Universal System units.
Definition: blobs.h:115
#define ZZ(v)
Definition: matrix1d.h:101
#define SPEED_UP_temps012
Definition: xmipp_macros.h:403

Variable Documentation

◆ blobs_conv_mutex

pthread_mutex_t blobs_conv_mutex = PTHREAD_MUTEX_INITIALIZER

Definition at line 31 of file blobs.cpp.

◆ slices_processed

int slices_processed

Definition at line 34 of file blobs.cpp.

◆ slices_status

int* slices_status

Definition at line 33 of file blobs.cpp.