458 const struct blobtype *blob = thread_data->blob;
461 int istep = thread_data->istep;
465 bool FORW = thread_data->FORW;
466 int eq_mode = thread_data->eq_mode;
468 int min_separation = thread_data->min_separation;
497 double vol_correction=0;
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) 511 bool condition = !FORW;
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;
529 for (
size_t i = 0; i < blob_table.
xdim; i++)
536 std::cout <<
"Blob (" << i <<
") r=" << (double)i / istep <<
537 " val= " <<
A1D_ELEM(blob_table, i) << std::endl;
553 return (
void*)
nullptr;
556 for(
int w = 0 ;
w < z_planes ;
w++ )
564 for(
int in = (
w-min_separation+1) ;
in <= (
w+min_separation-1 ) ;
in ++ )
568 if( (
in >= 0 ) && (
in < z_planes ))
581 while( assigned_slice == -1);
587 k = (int)(assigned_slice +
ZZ( grid->
lowest ));
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;
619 std::cout <<
"Center of the blob moved to " 622 << act_coord.transpose() << std::endl;
623 << real_position.transpose() << std::endl;
629 real_position = act_coord;
644 if (
XX(corner1) >=
xF)
646 if (
YY(corner1) >=
yF)
648 if (
ZZ(corner1) >=
zF)
650 if (
XX(corner2) <=
x0)
652 if (
YY(corner2) <=
y0)
654 if (
ZZ(corner2) <=
z0)
658 if (!process && condition)
659 std::cout <<
" It is outside output volume\n";
665 if (process && condition)
666 std::cout <<
" It is not interesting\n";
675 std::cout <<
"Corner 1 for this point " << corner1.transpose() << std::endl;
676 std::cout <<
"Corner 2 for this point " << corner2.transpose() << std::endl;
693 std::cout <<
"Clipped and rounded Corner 1 " << corner1.transpose() << std::endl;
694 std::cout <<
"Clipped and rounded Corner 2 " << corner2.transpose() << std::endl;
705 vol_correction = -1e38;
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++)
716 if (vol_mask !=
nullptr &&
A3D_ELEM(*vol_mask, intz, inty, intx)!=0.0)
720 VECTOR_R3(gcurrent, (
double)intx, (
double)inty, (
double)intz);
728 d =
sqrt(
XX(gcurrent) *
XX(gcurrent) +
729 YY(gcurrent) *
YY(gcurrent) +
730 ZZ(gcurrent) *
ZZ(gcurrent));
733 id = (int)(d * istep);
738 std::cout <<
"At (" << intx <<
"," 739 << inty <<
"," << intz <<
") distance=" <<
d;
748 A3D_ELEM(*vol_voxels, intz, inty, intx) +=
755 std::cout <<
" adding " <<
A3D_ELEM(*vol_blobs, k, i, j)
756 <<
" * " <<
A1D_ELEM(blob_table,
id) <<
" = " 758 A1D_ELEM(blob_table,
id) << std::endl;
762 if (vol_corr !=
nullptr)
763 A3D_ELEM(*vol_corr, intz, inty, intx) +=
768 double contrib =
A3D_ELEM(*vol_corr, intz, inty, intx) *
773 vol_correction += contrib;
777 if (contrib > vol_correction)
778 vol_correction = contrib;
785 std::cout <<
" adding " <<
A3D_ELEM(*vol_corr, intz, inty, intx)
786 <<
" * " <<
A1D_ELEM(blob_table,
id) <<
" = " 787 << contrib << std::endl;
798 A3D_ELEM(*vol_blobs, k, i, j) += vol_correction / N_eq;
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;
822 for(
int in = (assigned_slice-min_separation+1) ;
in <= (assigned_slice+min_separation-1);
in ++ )
824 if(
in != assigned_slice )
826 if( (
in >= 0 ) && (
in < z_planes))
838 return (
void*)
nullptr;
double kaiser_value(double r, double a, double alpha, int m)
double alpha
Smoothness parameter.
void resize(size_t Ndim, size_t Zdim, size_t Ydim, size_t Xdim, bool copy=true)
Matrix1D< double > highest
#define V3_PLUS_CT(a, b, c)
void grid2universe(const Matrix1D< double > &gv, Matrix1D< double > &uv) const
void sqrt(Image< double > &op)
#define V3_MINUS_V3(a, b, c)
pthread_mutex_t blobs_conv_mutex
void inv(Matrix2D< T > &result) const
bool is_interesting(const Matrix1D< double > &uv) const
Matrix1D< double > lowest
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
#define A3D_ELEM(V, k, i, j)
double relative_size
Measuring unit in the grid coordinate system.
#define M3x3_BY_V3x1(a, M, b)
#define VECTOR_R3(v, x, y, z)
int order
Derivation order and Bessel function order.
double radius
Spatial radius in Universal System units.
#define SPEED_UP_temps012