#include "DAGH.h"
/*
file Cluster3.C
date Sat Mar 23 17:37:08 1996
author Paul Walker
desc A 3D Clusterer using signature clustering.
There is also a 2-D interface
*/
/*
routine Cluster2
date Wed Apr 17 20:48:29 1996
author Paul Walker
desc A 2D interface to the recursive clustering algorithm.
*/
#include "DAGH.h"
void Cluster2(GridData(2)<short> &flag,
double Efficiency,
int MinWidth,
int BufferWidth,
BBoxList &Result) {
BBoxList Result3D;
BBox flat_3d(3,flag.bbox().lower(0),flag.bbox().lower(1),0,
flag.bbox().upper(0),flag.bbox().upper(1),0,
flag.bbox().stepsize(0),flag.bbox().stepsize(1),1);
GridData(2)<short> newflag(flag.bbox());
GridData(3)<short> clust_this(flat_3d);
void Cluster3R(GridData(3)<short> &flag, double, int,
BBoxList&);
int flagged_points = 0;
newflag = 0; clust_this = 0;
BBox buffer_bbox = flag.bbox();
buffer_bbox.grow(-BufferWidth);
for_2(i, j, buffer_bbox, flag.stepsize())
if (flag(i,j)) {
flagged_points ++;
BBox where(2,i,j,i,j,
flag.bbox().stepsize(0),
flag.bbox().stepsize(1));
where.grow(BufferWidth);
newflag.equals(1,where);
}
end_for;
if (flagged_points) {
for_2(i,j,flag.bbox(),flag.stepsize())
clust_this(i,j,0) = newflag(i,j);
end_for;
Cluster3R(clust_this, Efficiency, MinWidth, Result3D);
}
/*
Make sure that rank is appropriately set for the bboxes in
Result!!
Manish Parashar Thu May 9 07:54:32 CDT 1996
*/
for (BBox* tmpbb=Result3D.first();tmpbb;tmpbb=Result3D.next()) {
BBox tmp2d(2,tmpbb->lower(),tmpbb->upper(), tmpbb->stepsize());
Result.add(tmp2d);
}
}
/*
routine Cluster3
date Wed Apr 17 20:48:29 1996
author Paul Walker
desc
A 3D interface to the recursive clustering algorithm.
*/
void Cluster3(GridData(3)<short> &flag,
double Efficiency,
int MinWidth,
int BufferWidth,
BBoxList &Result) {
BBoxList tmpresult, result;
GridData(3)<short> newflag(flag.bbox());
void Cluster3R(GridData(3)<short> &flag, double, int,
BBoxList&);
int flagged_points = 0;
newflag = 0;
/* Buffer */
BBox buffer_bbox = flag.bbox();
buffer_bbox.grow(-BufferWidth);
for_3(i, j, k, buffer_bbox, flag.stepsize())
if (flag(i,j,k)) {
flagged_points ++;
BBox where(3,i,j,k,i,j,k,
flag.bbox().stepsize(0),
flag.bbox().stepsize(1),
flag.bbox().stepsize(2));
where.grow(BufferWidth);
newflag.equals(1,where);
}
end_for;
if (flagged_points) {
Cluster3R(newflag, Efficiency, MinWidth, Result);
}
}
/*
routine Cluster3R
date Sat Mar 23 18:35:36 1996
author Paul Walker
desc
A recursive signature clustering routine.
calls Cluster_Prune Cluster_Slice
*/
void Cluster3R(GridData(3)<short> &flag,
double Efficiency,
int MinWidth,
BBoxList &recurseOnThis) {
void Cluster3R(GridData(3)<short> &flag, double, int, BBoxList&);
void Cluster_Prune(GridData(1)<int> &sig, int *, int *, int);
void Cluster_Slice(GridData(1)<int> &sig, int *, int *, int);
int BiggestExtent = MinWidth;
BBoxList result, tmpresult;
int l[3], u[3], s[3];
int i;
for (i=0; i < 3; i++) s[i] = flag.bbox().stepsize(i);
// cout << "Clustering " << flag.bbox() <<"\n";
// cout.flush();
/* Declare and fill in Signature Planes in each direction */
Array(1)<int> i_sig(flag.bbox().lower(0),
flag.bbox().upper(0),
flag.bbox().stepsize(0));
Array(1)<int> j_sig(flag.bbox().lower(1),
flag.bbox().upper(1),
flag.bbox().stepsize(1));
Array(1)<int> k_sig(flag.bbox().lower(2),
flag.bbox().upper(2),
flag.bbox().stepsize(2));
i_sig = 0;
j_sig = 0;
k_sig = 0;
for_3(ii, jj, kk, flag.bbox(), flag.stepsize())
i_sig(ii) += flag(ii,jj,kk);
j_sig(jj) += flag(ii,jj,kk);
k_sig(kk) += flag(ii,jj,kk);
end_for;
/* Check for pruning */
int ipl, ipu, jpl, jpu, kpl, kpu;
Cluster_Prune(i_sig,&ipl,&ipu,BiggestExtent*flag.bbox().stepsize(0));
Cluster_Prune(j_sig,&jpl,&jpu,BiggestExtent*flag.bbox().stepsize(1));
Cluster_Prune(k_sig,&kpl,&kpu,BiggestExtent*flag.bbox().stepsize(2));
Coords pl(3,ipl,jpl,kpl);
Coords pu(3,ipu,jpu,kpu);
if (! ((pl == flag.bbox().lower()) &&
(pu == flag.bbox().upper()))) {
Coords ps = Coords(3,s);
BBox newB(pl,pu,ps);
GridData(3)<short> subgridflag(newB);
// cout <<" PRUNE :" << newB <<"\n";
// cout.flush();
subgridflag = flag;
/* Check my efficiency! */
int flagged = 0;
for_3(ii, jj, kk, subgridflag.bbox(), subgridflag.stepsize())
flagged += subgridflag(ii,jj,kk);
end_for
if (((double)flagged / (double)newB.size()) >= Efficiency) {
recurseOnThis.add(newB);
return;
} else {
Cluster3R(subgridflag,Efficiency, MinWidth, recurseOnThis);
return;
}
}
/* Find the slice positions */
int istr=-1, jstr=-1, kstr=-1;
int islice=-1, jslice=-1, kslice=-1;
Cluster_Slice(i_sig, &islice, &istr, BiggestExtent*flag.bbox().stepsize(0));
Cluster_Slice(j_sig, &jslice, &jstr, BiggestExtent*flag.bbox().stepsize(1));
Cluster_Slice(k_sig, &kslice, &kstr, BiggestExtent*flag.bbox().stepsize(2));
// cout << "SLICE:" << islice << " " <<
// jslice << " " << kslice << "\n";
// cout << "STR :" << istr << " " <<
// jstr << " " << kstr << "\n";
// cout.flush();
if (islice+jslice+kslice == -3) {
// cout << "RETURNING " << flag.bbox() << endl
// << flush;
recurseOnThis.add(flag.bbox());
return;
}
/* Split along the line with the biggest strength */
BBox divme = flag.bbox();
if (istr > jstr && istr > kstr) {
l[0] = divme.lower(0); l[1]= divme.lower(1); l[2] = divme.lower(2);
u[0] = islice; u[1]= divme.upper(1); u[2] = divme.upper(2);
tmpresult.add(BBox(3,l,u,s));
l[0] = islice+s[0]; l[1]= divme.lower(1); l[2] = divme.lower(2);
u[0] = divme.upper(0); u[1]= divme.upper(1); u[2] = divme.upper(2);
tmpresult.add(BBox(3,l,u,s));
} else if (jstr > kstr) {
l[0] = divme.lower(0); l[1]= divme.lower(1); l[2] = divme.lower(2);
u[0] = divme.upper(0); u[1]= jslice; u[2] = divme.upper(2);
tmpresult.add(BBox(3,l,u,s));
l[0] = divme.lower(0); l[1]= jslice+s[1]; l[2] = divme.lower(2);
u[0] = divme.upper(0); u[1]= divme.upper(1); u[2] = divme.upper(2);
tmpresult.add(BBox(3,l,u,s));
} else {
l[0] = divme.lower(0); l[1]= divme.lower(1); l[2] = divme.lower(2);
u[0] = divme.upper(0); u[1]= divme.upper(1); u[2] = kslice;
tmpresult.add(BBox(3,l,u,s));
l[0] = divme.lower(0); l[1]= divme.lower(1); l[2] = kslice+s[2];
u[0] = divme.upper(0); u[1]= divme.upper(1); u[2] = divme.upper(2);
tmpresult.add(BBox(3,l,u,s));
}
/* foreach New bbox */
for (BBox* newB = tmpresult.first(); newB; newB = tmpresult.next()) {
/* Figure out how many points are flagged */
GridData(3)<short> subgridflag(*newB);
// cout << "Handling " << *newB << "\n";
subgridflag = flag;
int flagged = 0;
double myEff;
for_3(ii, jj, kk, subgridflag.bbox(), subgridflag.stepsize())
// cout << "SGF: "<< ii <<" "<< jj <<" "
// << kk <<" " <<subgridflag(ii,jj,kk)<<"\n";
// cout.flush();
flagged += subgridflag(ii,jj,kk);
end_for
myEff = (double)flagged/(double)((*newB).size());
/* Do we have ANY? */
if (flagged) {
// cout <<"Flagged " << flagged << " "
// << *newB << "\n";
// cout <<"Efficiency "<< myEff << " of "
// << (*newB).size() <<"\n";
// cout.flush();
/* Stopping Condition */
if (myEff >= Efficiency ||
(((*newB).size()) < 8) ||
(islice + jslice + kslice == -3)) {
/* No. This one stays = the list if it contains flagged points */
/* Don't add a degenerate grid! */
int add = 1;
for (BBox *tmp = recurseOnThis.first(); tmp;
tmp = recurseOnThis.next()) {
if ((*tmp).inside((*newB).lower()) &&
(*tmp).inside((*newB).upper())) add = 0;
}
if (add) {
// cout << " Adding "<< *newB <<"\n";
// cout.flush();
recurseOnThis.add(*newB);
}
} else {
/* Yes: Recluster this sub-box */
// cout << "FALL: Flagged " << flagged << " "
// << *newB << "\n";
// cout <<"FALL: Efficiency "<<myEff <<" of "
// << (*newB).size() <<"\n";
// cout.flush();
BBoxList subcluster;
Cluster3R(subgridflag, Efficiency,
MinWidth, subcluster);
for (BBox *me = subcluster.first(); me; me = subcluster.next())
recurseOnThis.add(*me);
}
}
}
return;
}
/*
routine Cluster_Prune
date Mon Mar 25 15:12:06 1996
author Paul Walker
desc
Prunes out external zeros and returns a lower and upper bound for
prune locations.
*/
void Cluster_Prune(GridData(1)<int> &sig, int *pl, int *pu,
int bw) {
*pl = sig.bbox().lower(0);
*pu = sig.bbox().upper(0);
if (*pu - *pl <= bw) {
return;
}
while (sig(*pl) == 0 && *pl < sig.bbox().upper(0))
*pl += sig.bbox().stepsize(0);
while (sig(*pu) == 0 && *pu > sig.bbox().lower(0))
*pu -= sig.bbox().stepsize(0);
if (*pu - *pl < bw) {
int diff = bw - (*pu - *pl);
int addl, addu;
int udist = sig.bbox().upper(0) - *pu;
int ldist = *pl - sig.bbox().lower(0);
addl = ((diff/sig.bbox().stepsize(0))/2) * sig.bbox().stepsize(0);
addu = diff - addl;
if (addu > udist) {
addu = udist;
addl = diff - addu;
}
if (addl > ldist) {
addl = ldist;
addu = diff - addu;
}
*pu += addu;
*pl -= addl;
}
}
/*
routine Cluster_Slice
date Mon Mar 25 15:09:21 1996
author Paul Walker
desc
Looks for a laplacian or zero slice in the internal grid area.
Since prune has already been called, we know this will happen
inside!
*/
void Cluster_Slice(GridData(1)<int> &sig, int *slice, int *str,
int bw) {
int i, l, u, s, t, a; // ctr, low, up, step, str, answer
l = sig.bbox().lower(0);
u = sig.bbox().upper(0);
s = sig.bbox().stepsize(0);
GridData(1)<int> lap(l,u,s);
/* Since this box is already pruned, we don't want to split it if
it is less than 2 the min width */
if (u-l < 2*bw) {
*str = -1;
*slice = -1;
return;
}
t = -1;
a = -1;
/* Look for a zero near the center */
for (i=l; i <= u-s; i+=s) {
if (sig(i) == 0) {
int di = (i-l < u-i ? i-l : u-i);
if (di > t) {
a = i;
t = di;
}
}
}
if (a != -1) {
*slice = a;
*str = 10000;
return;
}
t = -1;
a = -1;
/* Evaluate the laplacian etc */
double malap = 0;
for (i = l+s; i <= u-s; i+=s) {
lap(i) = sig(i+s) + sig(i-s) - 2*sig(i);
if (abs(lap(i)) > malap) malap = abs(lap(i));
}
/* And find the strongest second derivative */
if (malap > 0) {
// If it isn't then lap is zero everywhere so don't slice either...
// This case isn't included here because of the (needed) <=
// for the graze-zero-but-interesting case which still makes it
for (i=l+bw; i<=u-bw; i+=s)
if (lap(i+s)*lap(i) <= 0) {
int ts = abs(lap(i+s)-lap(i));
if (ts > t) {
t = ts;
a = i;
// cout << "PICKING slice " << i << endl
// << flush;
// cout << ts << " " << t << " "
// << lap(i+s) << " " << lap(i) <<
// endl << flush;
}
}
}
*slice = a;
*str = t;
}