Skip to content

Commit

Permalink
move Graphics::TriD::Rout::{combcoords,attract,repulse} -> ImageND, test
Browse files Browse the repository at this point in the history
  • Loading branch information
@mohawk2
mohawk2 committed Dec 3, 2024
1 parent 4d28650 commit 5901118
Show file tree
Hide file tree
Showing 6 changed files with 192 additions and 192 deletions.
160 changes: 160 additions & 0 deletions Basic/lib/PDL/ImageND.pd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1192,4 +1192,164 @@ sub PDL::path_segs {
EOPM
pp_add_exported('','path_segs');

pp_def('combcoords',
GenericTypes => ['F','D'],
DefaultFlow => 1,
Pars => 'x(); y(); z();
float [o]coords(tri=3);',
Code => '
$coords(tri => 0) = $x();
$coords(tri => 1) = $y();
$coords(tri => 2) = $z();
',
Doc => <<EOT
=for ref
Combine three coordinates into a single ndarray.
Combine x, y and z to a single ndarray the first dimension
of which is 3. This routine does dataflow automatically.
EOT
);

pp_def('repulse',
GenericTypes => ['F','D'],
Pars => 'coords(nc,np); [o]vecs(nc,np); int [t]links(np);',
OtherPars => '
double boxsize;
int dmult;
double a;
double b;
double c;
double d;
',
Code => <<'EOF',
double a = $COMP(a);
double b = $COMP(b);
double c = $COMP(c);
double d = $COMP(d);
int ind; int x,y,z;
HV *hv = newHV();
double boxsize = $COMP(boxsize);
int dmult = $COMP(dmult);
loop(np) %{
int index = 0;
$links() = -1;
loop(nc) %{
$vecs() = 0;
index *= dmult;
index += (int)($coords()/boxsize);
%}
/* Repulse old (shame to use x,y,z...) */
for (x=-1; x<=1; x++) {
for (y=-1; y<=1; y++) {
for (z=-1; z<=1; z++) {
int ni = index + x + dmult * y + dmult * dmult * z;
SV **svp = hv_fetch(hv, (char *)&ni, sizeof(int), 0);
if (svp && *svp) {
ind = SvIV(*svp) - 1;
while (ind>=0) {
double dist = 0;
double dist2;
double tmp;
double func;
loop(nc) %{
tmp = ($coords() - $coords(np => ind));
dist += tmp * tmp;
%}
dist = sqrt(1/(sqrt(dist)+d));
func = c * dist;
dist2 = dist * dist;
func += b * dist2;
dist2 *= dist2;
func += a * dist2;
loop(nc) %{
tmp = ($coords() - $coords(np => ind));
$vecs() -= func * tmp;
$vecs(np => ind) += func * tmp;
%}
ind = $links(np => ind);
}
}
}
}
}
/* Store new */
SV **svp = hv_fetch(hv, (char *)&index, sizeof(index), 1);
if(!svp || !*svp)
$CROAK("Invalid sv from hvfetch");
SV *sv = *svp;
int npv;
if(SvOK(sv) && (npv = SvIV(sv))) {
npv --;
$links() = $links(np => npv);
$links(np => npv) = np;
} else {
sv_setiv(sv,np+1);
$links() = -1;
}
%}
hv_undef(hv);
EOF
Doc => <<'EOF',
=for ref
Repulsive potential for molecule-like constructs.
C<repulse> uses a hash table of cubes to quickly calculate
a repulsive force that vanishes at infinity for many
objects. For use by the module L<PDL::Graphics::TriD::MathGraph>.
Checks all neighbouring boxes. The formula is:
(r = |dist|+d) a*r^-2 + b*r^-1 + c*r^-0.5
EOF
);

pp_def('attract',
GenericTypes => ['F','D'],
Pars => 'coords(nc,np);
int from(nl);
int to(nl);
strength(nl);
[o]vecs(nc,np);',
OtherPars => '
double m;
double ms;
',
Code => <<'EOF',
double m = $COMP(m);
double ms = $COMP(ms);
loop(nc,np) %{ $vecs() = 0; %}
loop(nl) %{
int f = $from();
int t = $to();
double s = $strength();
double dist = 0;
double tmp;
loop(nc) %{
tmp = $coords(np => f) - $coords(np => t);
dist += tmp * tmp;
%}
s *= ms * dist + m * sqrt(dist);
loop(nc) %{
tmp = $coords(np => f) - $coords(np => t);
$vecs(np => f) -= tmp * s;
$vecs(np => t) += tmp * s;
%}
%}
EOF
Doc => '
=for ref
Attractive potential for molecule-like constructs.
C<attract> is used to calculate
an attractive force for many
objects, of which some attract each other (in a way
like molecular bonds).
For use by the module L<PDL::Graphics::TriD::MathGraph>.
For definition of the potential, see the actual function.
'
);

pp_done();
22 changes: 22 additions & 0 deletions Basic/t/imagend.t
View file
Original file line number Diff line number Diff line change
Expand Up @@ -124,4 +124,26 @@ for (
is_deeply [map "$_", path_segs($pi, $p)], ['[0 1 2 3 1]', '[2 3]', '[4 5]'];
}

{
my ($x, $y, $z) = map float($_), 5..7;
my $c3 = combcoords($x,$y,$z);
is_pdl $c3, float(5,6,7);
$x++;
is_pdl $c3, float(6,6,7);
}

{
my $coords = float([0,-1,0], [-1,-1,-2], [3,5,2], [2,1,-3], [1,3,1], [1,1,2]);
my $from = indx([0,1,2,3,4,4,4,5,5,5]);
my $to = indx([1,2,3,1,0,2,3,0,1,2]);
is_pdl repulse($coords,3,5000,-100,-5,-0.1,0.01),
float('[-19.9785 -83.469 27.7307; -75.2499 -57.671 -78.2403;
57.7644 92.1679 42.674; 47.2407 6.71818 -93.5201;
-15.8246 75.1071 8.60458; 6.04795 -32.8531 92.7511]'), {atol=>1e-2};
is_pdl attract($coords,$from,$to,1,30,1),
float('[172.197 779.117 199.115; 2054.32 2486.76 1963.34;
-2004.3 -3652.66 -2542.66; -300.804 1010.14 1942.27;
211.198 -700.071 -680.188; -132.611 76.7175 -881.878]'), {atol=>1e-2};
}

done_testing;
1 change: 1 addition & 0 deletions Changes
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,6 +19,7 @@
- split PDL::Perldl2 out to separate distro
- fixed overly-rigid parsing of ASCII STL files (#504) - thanks Shugo for report
- split PDL::Complex out to separate distro
- move combcoords attract repulse from Graphics::TriD::Rout to ImageND

2.095 2024-11-03
- add PDL_GENTYPE_IS_{REAL,FLOATREAL,COMPLEX,SIGNED,UNSIGNED}_##ppsym (#502)
Expand Down
176 changes: 1 addition & 175 deletions Graphics/TriD/Rout/rout.pd
View file
Original file line number Diff line number Diff line change
Expand Up @@ -12,184 +12,10 @@ PDL::Graphics::TriD::Rout - Helper routines for Three-dimensional graphics
=head1 DESCRIPTION
This module is for miscellaneous PP-defined utility routines for
the PDL::Graphics::TriD module. Currently, there are
the PDL::Graphics::TriD module.
EOD

pp_def(
'combcoords',
GenericTypes => ['F','D'],
DefaultFlow => 1,
Pars => 'x(); y(); z();
float [o]coords(tri=3);',
Code => '
$coords(tri => 0) = $x();
$coords(tri => 1) = $y();
$coords(tri => 2) = $z();
',
Doc => <<EOT
=for ref
Combine three coordinates into a single ndarray.
Combine x, y and z to a single ndarray the first dimension
of which is 3. This routine does dataflow automatically.
EOT

);

# checks all neighbouring boxes.
# Returns (r = |dist|+d) a*r^-2 + b*r^-1 + c*r^-0.5
pp_def(
'repulse',
GenericTypes => ['F','D'],
Pars => 'coords(nc,np);
[o]vecs(nc,np);
int [t]links(np);',
OtherPars => '
double boxsize;
int dmult;
double a;
double b;
double c;
double d;
',
Code => '
double a = $COMP(a);
double b = $COMP(b);
double c = $COMP(c);
double d = $COMP(d);
int ind; int x,y,z;
HV *hv = newHV();
double boxsize = $COMP(boxsize);
int dmult = $COMP(dmult);
loop(np) %{
int index = 0;
$links() = -1;
loop(nc) %{
$vecs() = 0;
index *= dmult;
index += (int)($coords()/boxsize);
%}
/* Repulse old (shame to use x,y,z...) */
for(x=-1; x<=1; x++) {
for(y=-1; y<=1; y++) {
for(z=-1; z<=1; z++) {
int ni = index + x + dmult * y +
dmult * dmult * z;
SV **svp = hv_fetch(hv, (char *)&ni, sizeof(int),
0);
if(svp && *svp) {
ind = SvIV(*svp) - 1;
while(ind>=0) {
double dist = 0;
double dist2;
double tmp;
double func;
loop(nc) %{
tmp =
($coords() -
$coords(np => ind));
dist += tmp * tmp;
%}
dist = sqrt(1/(sqrt(dist)+d));
func = c * dist;
dist2 = dist * dist;
func += b * dist2;
dist2 *= dist2;
func += a * dist2;
loop(nc) %{
tmp =
($coords() -
$coords(np => ind));
$vecs() -=
func * tmp;
$vecs(np => ind) +=
func * tmp;
%}
ind = $links(np => ind);
}
}
}
}
}
/* Store new */
SV **svp = hv_fetch(hv, (char *)&index, sizeof(index), 1);
if(!svp || !*svp)
$CROAK("Invalid sv from hvfetch");
SV *sv = *svp;
int npv;
if(SvOK(sv) && (npv = SvIV(sv))) {
npv --;
$links() = $links(np => npv);
$links(np => npv) = np;
} else {
sv_setiv(sv,np+1);
$links() = -1;
}
%}
hv_undef(hv);
', Doc => '
=for ref
Repulsive potential for molecule-like constructs.
C<repulse> uses a hash table of cubes to quickly calculate
a repulsive force that vanishes at infinity for many
objects. For use by the module L<PDL::Graphics::TriD::MathGraph>.
For definition of the potential, see the actual function.
'
);

pp_def(
'attract',
GenericTypes => ['F','D'],
Pars => 'coords(nc,np);
int from(nl);
int to(nl);
strength(nl);
[o]vecs(nc,np);',
OtherPars => '
double m;
double ms;
',
Code => '
double m = $COMP(m);
double ms = $COMP(ms);
loop(nc,np) %{ $vecs() = 0; %}
loop(nl) %{
int f = $from();
int t = $to();
double s = $strength();
double dist = 0;
double tmp;
loop(nc) %{
tmp = $coords(np => f) -
$coords(np => t);
dist += tmp * tmp;
%}
s *= ms * dist + m * sqrt(dist);
loop(nc) %{
tmp = $coords(np => f) -
$coords(np => t);
$vecs(np => f) -= tmp * s;
$vecs(np => t) += tmp * s;
%}
%}
', Doc => '
=for ref
Attractive potential for molecule-like constructs.
C<attract> is used to calculate
an attractive force for many
objects, of which some attract each other (in a way
like molecular bonds).
For use by the module L<PDL::Graphics::TriD::MathGraph>.
For definition of the potential, see the actual function.
'
);

sub trid {
my ($par,$ind) = @_;
join ',', map {"\$$par($ind => $_)"} (0..2);
Expand Down
Loading

0 comments on commit 5901118

Please sign in to comment.