package edu.stanford.rsl.conrad.volume3d;
import java.util.Arrays;
public class MaxEigenValue implements Runnable {
private int x = -1, y, z;
private float [][][] T;
private int dimensions;
private float [] eigmax;
private boolean finished = true;
private boolean data = false;
private boolean terminate = false;
public MaxEigenValue(int dimensions){
this.dimensions = dimensions;
data = false;
finished = true;
}
public void terminate (){
terminate = true;
}
public void setData(int x, int y, int z, float [][][] T){
this.x = x;
this.y = y;
this.z = z;
this.T = T;
data = true;
finished = false;
eigmax = new float[z];
}
public boolean done(){
return finished;
}
public void getData(float [][][] array){
if (x >= 0){
for (int i = 0; i< z; i++){
array[x][y][i] = eigmax[i];
}
}
}
private void nrerror(Object s){
//TODO: nothing
//System.out.println("Jacobi " +x + " " + y + " "+ z + ": " +s.toString());
}
private static void ROTATE(float [][] a, int i, int j, int k, int l, float s, float tau) {
float g = a[i][j];
float h = a[k][l];
a[i][j]=g-s*(h+g*tau);
a[k][l]=h+s*(g-h*tau);
}
private static float fabs(float in){
return Math.abs(in);
}
private static float sqrt(double in){
return (float) Math.sqrt(in);
}
public static float [] vector(int i, int n){
float [] revan = new float [n];
Arrays.fill(revan, i);
return revan;
}
public static void jacobi(float [][] a, int n, float []d, float [][]v, Integer nrot){
int j,iq,ip,i;
float tresh,theta,tau,t,sm,s,h,g,c;
float [] b;
float [] z;
b=vector(1,n);
z=vector(1,n);
for (ip=0;ip<n;ip++) {
for (iq=0;iq<n;iq++) v[ip][iq]=0.0f;
v[ip][ip]=1.0f;
}
for (ip=0;ip<n;ip++) {
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0f;
}
nrot=0;
for (i=0;i<50;i++) {
sm=0.0f;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++)
sm += Math.abs(a[ip][iq]);
}
if (sm == 0.0) {
z=null;
b=null;
return;
}
if (i < 4)
tresh= (0.2f*sm/(n*n));
else
tresh=0.0f;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++) {
g=100.0f*Math.abs(a[ip][iq]);
if (i > 4 && fabs(d[ip])+g == fabs(d[ip])
&& fabs(d[iq])+g == fabs(d[iq]))
a[ip][iq]=0.0f;
else if (fabs(a[ip][iq]) > tresh) {
h=d[iq]-d[ip];
if (fabs(h)+g == fabs(h))
t=(a[ip][iq])/h;
else {
theta=0.5f*h/(a[ip][iq]);
t=1.0f/(fabs(theta)+sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0f/sqrt(1+t*t);
s=t*c;
tau=s/(1.0f+c);
h=t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip][iq]=0.0f;
for (j=0;j<=ip-1;j++) {
ROTATE(a,j,ip,j,iq,s,tau);
g = a[j][ip];
h = a[j][iq];
}
for (j=ip+1;j<=iq-1;j++) {
ROTATE(a,ip,j,j,iq,s,tau);
g = a[ip][j];
h = a[j][iq];
}
for (j=iq+1;j<n;j++) {
ROTATE(a,ip,j,iq,j,s,tau);
g = a[ip][j];
h = a[iq][j];
}
for (j=0;j<n;j++) {
ROTATE(v,j,ip,j,iq,s,tau);
g = v[j][ip];
h = v[j][iq];
}
++(nrot);
}
}
}
for (ip=0;ip<n;ip++) {
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0f;
}
}
//nrerror("Too many iterations in routine JACOBI");
}
public void run() {
while (!terminate){
if(data){
for (int i = 0; i < z; i++){
float [][] v = new float[Volume3D.MAX_DIM][Volume3D.MAX_DIM];
float [] eig = new float [Volume3D.MAX_DIM];
int dim_loop;
Integer nrot = new Integer(0);
// TODO Auto-generated method stub
jacobi(T[i], dimensions, eig, v, nrot);
eigmax[i]=eig[0];
for (dim_loop=1; dim_loop<dimensions; dim_loop++)
if (eig[dim_loop]>eigmax[i]) eigmax[i]=eig[dim_loop];
nrerror("Maximal Value:" + eigmax[i]);
}
finished = true;
data = false;
} else {
try {
Thread.sleep(0);
} catch (InterruptedException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
}
}
}
/*
* Copyright (C) 2010-2014 Andreas Maier
* CONRAD is developed as an Open Source project under the GNU General Public License (GPL).
*/