Add vector acceleration on Lucy and Lucy-TV !

.classpath

@@ -6,7 +6,7 @@
 			<attribute name="maven.pomderived" value="true"/>
 		</attributes>
 	</classpathentry>
-	<classpathentry excluding="**|ij/" kind="src" output="target/classes" path="src/main/resources">
+	<classpathentry excluding="**" kind="src" output="target/classes" path="src/main/resources">
 		<attributes>
 			<attribute name="maven.pomderived" value="true"/>
 		</attributes>
@@ -23,5 +23,20 @@
 	</classpathentry>
 	<classpathentry kind="lib" path="icy.jar"/>
 	<classpathentry kind="lib" path="jtransforms.jar"/>
+	<classpathentry kind="src" output="target/test-classes" path="src/test/java">
+		<attributes>
+			<attribute name="optional" value="true"/>
+			<attribute name="maven.pomderived" value="true"/>
+		</attributes>
+	</classpathentry>
+	<classpathentry kind="lib" path="C:/Users/Administrator/.m2/repository/edu/emory/mathcs/jtransforms/2.4/jtransforms-2.4.jar">
+		<attributes>
+			<attribute name="maven.pomderived" value="true"/>
+			<attribute name="maven.groupId" value="edu.emory.mathcs"/>
+			<attribute name="maven.artifactId" value="jtransforms"/>
+			<attribute name="maven.version" value="2.4"/>
+			<attribute name="maven.scope" value="runtime"/>
+		</attributes>
+	</classpathentry>
 	<classpathentry kind="output" path="target/classes"/>
 </classpath>

.settings/org.eclipse.core.resources.prefs

@@ -1,2 +1,4 @@
 eclipse.preferences.version=1
+encoding//src/main/java=UTF-8
+encoding//src/main/resources=UTF-8
 encoding/<project>=UTF-8

src/main/java/bilib/component/PanelImage.java

@@ -33,6 +33,7 @@
 import java.awt.Graphics;
 import java.awt.Image;
 import java.awt.image.BufferedImage;
+
 import javax.swing.JPanel;
 
 import resources.ImageLoader;

src/main/java/deconvolution/algorithm/RichardsonLucy.java

@@ -50,23 +50,108 @@ public RichardsonLucy(int iterMax) {
 	public RealSignal call() {
 		ComplexSignal H = fft.transform(h);
 		ComplexSignal U = new ComplexSignal("RL-U", y.nx, y.ny, y.nz);
-		RealSignal x = y.duplicate();
+		RealSignal x = new RealSignal("x",y.nx,y.ny,y.nz);
 		RealSignal p = y.duplicate();
 		RealSignal u = y.duplicate();
+
+		// For vector acceleration	
+		RealSignal yv_update = new RealSignal("yv_update",y.nx,y.ny,y.nz);
+		RealSignal v_vector = new RealSignal("v_vector",y.nx,y.ny,y.nz);	
+		RealSignal v_vector1 = new RealSignal("v_vector",y.nx,y.ny,y.nz);	
+		RealSignal xv_update = new RealSignal("xv_update",y.nx,y.ny,y.nz);
+		RealSignal y_vector = y.duplicate();
+		RealSignal vv_update = v_vector.duplicate();
+		RealSignal x_update = x.duplicate();
+		//First iteration ! 
+		float alpha1 = 0;
+		float alphal1 = 0;
+		float alphau1 = 0;		
+		fft.transform(y_vector, U);
+		U.times(H);
+		fft.inverse(U, u);
+		Operations.divide(y, u, p);
+		fft.transform(p, U);
+		U.timesConjugate(H);
+		fft.inverse(U, u);
+		for (int z = 0; z < y.nz; z++) {
+			for (int i = 0; i < y.nx * y.ny; i++) {
+				System.out.print(i);
+				x.data[z][i]= y_vector.data[z][i] * u.data[z][i];
+			}
+		}
+		Operations.subtract(x, y_vector, v_vector);	
+		//prepare for vector
+		fft.transform(yv_update, U);
+		U.times(H);
+		fft.inverse(U, u);
+		Operations.divide(y, u, p);
+		fft.transform(p, U);
+		U.timesConjugate(H);
+		fft.inverse(U, u);
+		for (int z = 0; z < y.nz; z++) {
+			for (int i = 0; i < y.nx * y.ny; i++) {
+				xv_update.data[z][i]= yv_update.data[z][i] * u.data[z][i];
+			}
+		}
+		Operations.subtract(xv_update, yv_update, v_vector1);
+		for (int z = 0; z < y.nz; z++) {
+			for (int i = 0; i < y.nx * y.ny; i++) {
+				alphau1+= vv_update.data[z][i] * v_vector1.data[z][i];
+				alphal1 += v_vector1.data[z][i] * v_vector1.data[z][i]+0.001;
+			}
+		}
+		alpha1=alphau1/alphal1;
+		if (alpha1<0)
+			alpha1=(float) 0.0001;
+		if (alpha1>1)
+			alpha1=1;
+		y_vector=Operations.subtract(x, x_update).times(alpha1).plus(x);		
+
+		//iteration really start!
 		while (!controller.ends(x)) {
-			fft.transform(x, U);
+			float alpha = 0;
+			float alphal = 0;
+			float alphau = 0;
+			// For vector acceleration
+			x_update = x.duplicate();
+			fft.transform(y_vector, U);
 			U.times(H);
 			fft.inverse(U, u);
 			Operations.divide(y, u, p);
 			fft.transform(p, U);
 			U.timesConjugate(H);
 			fft.inverse(U, u);
-			x.times(u);
+			for (int z = 0; z < y.nz; z++) {
+				for (int i = 0; i < y.nx * y.ny; i++) {
+					x.data[z][i]= y_vector.data[z][i] * u.data[z][i];
+				}
+			}
+			vv_update = v_vector.duplicate();
+			Operations.subtract(x, y_vector, v_vector);
+			for (int z = 0; z < y.nz; z++) {
+				for (int i = 0; i < y.nx * y.ny; i++) {
+					alphau+= vv_update.data[z][i] * v_vector.data[z][i];
+					alphal += vv_update.data[z][i] * vv_update.data[z][i]+0.001;
+				}
+			}
+			alpha=alphau/alphal;
+			if (alpha<0)
+				alpha=(float) 0.0001;
+			if (alpha>1)
+				alpha=1;
+			System.out.print(alpha);
+			y_vector=Operations.subtract(x, x_update).times(alpha).plus(x);
 		}
 		SignalCollector.free(H);
 		SignalCollector.free(p);
 		SignalCollector.free(u);
 		SignalCollector.free(U);
+		SignalCollector.free(yv_update);
+		SignalCollector.free(v_vector );
+		SignalCollector.free(v_vector1);
+		SignalCollector.free(xv_update);
+		SignalCollector.free(y_vector);
+		SignalCollector.free(vv_update);
 		return x;
 	}
 

src/main/java/deconvolution/algorithm/RichardsonLucyTV.java

@@ -60,21 +60,96 @@ public RealSignal call() {
 		RealSignal ggx = y.duplicate();
 		RealSignal ggy = y.duplicate();
 		RealSignal ggz = y.duplicate();
-
 		RealSignal u  = gx;	// resued memory
 		RealSignal p  = gy;	// resued memory
 		RealSignal tv = gz; // resued memory
+
+		// For vector acceleration	
+		RealSignal yv_update = new RealSignal("yv_update",y.nx,y.ny,y.nz);
+		RealSignal v_vector = new RealSignal("v_vector",y.nx,y.ny,y.nz);	
+		RealSignal v_vector1 = new RealSignal("v_vector",y.nx,y.ny,y.nz);	
+		RealSignal xv_update = new RealSignal("xv_update",y.nx,y.ny,y.nz);
+		RealSignal y_vector = y.duplicate();
+		RealSignal vv_update = v_vector.duplicate();
+		RealSignal x_update = x.duplicate();
+		//First iteration ! 
+		float alpha1 = 0;
+		float alphal1 = 0;
+		float alphau1 = 0;
+		gradientX(y_vector, gx);
+		gradientY(y_vector, gy);
+		gradientZ(y_vector, gz);
+		normalize(gx, gy, gz);	
+		gradientX(gx, ggx);
+		gradientY(gy, ggy);
+		gradientZ(gz, ggz);
+		compute((float)lambda, ggx, ggy, ggz, tv);
+		fft.transform(y_vector, U);
+		U.times(H);
+		fft.inverse(U, u);
+		Operations.divide(y, u, p);
+		fft.transform(p, U);
+		U.timesConjugate(H);
+		fft.inverse(U, u);
+		for (int z = 0; z < y.nz; z++) {
+			for (int i = 0; i < y.nx * y.ny; i++) {
+				System.out.print(i);
+				x.data[z][i]= y_vector.data[z][i] * u.data[z][i]*tv.data[z][i];
+			}
+		}
+		Operations.subtract(x, y_vector, v_vector);	
+		//prepare for vector
+		gradientX(yv_update, gx);
+		gradientY(yv_update, gy);
+		gradientZ(yv_update, gz);
+		normalize(gx, gy, gz);	
+		gradientX(gx, ggx);
+		gradientY(gy, ggy);
+		gradientZ(gz, ggz);
+		compute((float)lambda, ggx, ggy, ggz, tv);
+		fft.transform(yv_update, U);
+		U.times(H);
+		fft.inverse(U, u);
+		Operations.divide(y, u, p);
+		fft.transform(p, U);
+		U.timesConjugate(H);
+		fft.inverse(U, u);
+		for (int z = 0; z < y.nz; z++) {
+			for (int i = 0; i < y.nx * y.ny; i++) {
+				xv_update.data[z][i]= yv_update.data[z][i] * u.data[z][i]*tv.data[z][i];
+			}
+		}
+		Operations.subtract(xv_update, yv_update, v_vector1);
+		for (int z = 0; z < y.nz; z++) {
+			for (int i = 0; i < y.nx * y.ny; i++) {
+				alphau1+= vv_update.data[z][i] * v_vector1.data[z][i];
+				alphal1 += v_vector1.data[z][i] * v_vector1.data[z][i]+0.001;
+			}
+		}
+		alpha1=alphau1/alphal1;
+		if (alpha1<0)
+			alpha1=(float) 0.0001;
+		if (alpha1>1)
+			alpha1=1;
+		y_vector=Operations.subtract(x, x_update).times(alpha1).plus(x);		
 
+
+
 		while(!controller.ends(x)) {
-			gradientX(x, gx);
-			gradientY(x, gy);
-			gradientZ(x, gz);
+			float alpha = 0;
+			float alphal = 0;
+			float alphau = 0;
+			// For vector acceleration
+			x_update = x.duplicate();	
+			gradientX(y_vector, gx);
+			gradientY(y_vector, gy);
+			gradientZ(y_vector, gz);
 			normalize(gx, gy, gz);	
 			gradientX(gx, ggx);
 			gradientY(gy, ggy);
 			gradientZ(gz, ggz);
 			compute((float)lambda, ggx, ggy, ggz, tv);
-			fft.transform(x, U);
+			fft.transform(y_vector, U);
 			U.times(H);
 			fft.inverse(U, u);
 			Operations.divide(y, u, p);
@@ -83,6 +158,21 @@ public RealSignal call() {
 			fft.inverse(U, u);
 			x.times(u); 
 			x.times(tv);
+			// For vector acceleration
+			vv_update = v_vector.duplicate();
+			Operations.subtract(x, y_vector, v_vector);
+			for (int z = 0; z < y.nz; z++) {
+				for (int i = 0; i < y.nx * y.ny; i++) {
+					alphau += vv_update.data[z][i] * v_vector.data[z][i];
+					alphal += vv_update.data[z][i] * vv_update.data[z][i];
+				}
+			}
+			alpha=alphau/alphal;
+			if (alpha<0)
+				alpha=(float) 0.0001;
+			if (alpha>1)
+				alpha=1;
+			y_vector=Operations.subtract(x, x_update).times(alpha).plus(x);	
 		}
 		SignalCollector.free(H);
 		SignalCollector.free(U);
@@ -92,6 +182,12 @@ public RealSignal call() {
 		SignalCollector.free(tv);
 		SignalCollector.free(u);
 		SignalCollector.free(p);
+		SignalCollector.free(yv_update);
+		SignalCollector.free(v_vector );
+		SignalCollector.free(v_vector1);
+		SignalCollector.free(xv_update);
+		SignalCollector.free(y_vector);
+		SignalCollector.free(vv_update);
 		return x;
 	}
 

src/main/java/deconvolution/algorithm/Simulation.java

@@ -33,11 +33,11 @@
 
 import java.util.concurrent.Callable;
 
+import bilib.tools.PsRandom;
 import signal.ComplexSignal;
 import signal.Operations;
 import signal.RealSignal;
 import signal.SignalCollector;
-import bilib.tools.PsRandom;
 
 public class Simulation extends Algorithm implements Callable<RealSignal> {
 

src/main/java/deconvolution/capsule/AlgorithmCapsule.java

@@ -37,9 +37,6 @@
 import javax.swing.JSplitPane;
 import javax.swing.JViewport;
 
-import signal.ComplexSignal;
-import signal.RealSignal;
-import signal.SignalCollector;
 import bilib.component.HTMLPane;
 import bilib.table.CustomizedTable;
 import bilib.tools.NumFormat;
@@ -51,6 +48,9 @@
 import deconvolutionlab.monitor.Monitors;
 import fft.AbstractFFT;
 import fft.FFT;
+import signal.ComplexSignal;
+import signal.RealSignal;
+import signal.SignalCollector;
 
 /**
  * This class is a information module for the algorithm.

src/main/java/deconvolution/capsule/ImageCapsule.java

@@ -31,11 +31,8 @@
 
 package deconvolution.capsule;
 
-import java.awt.Dimension;
-
 import javax.swing.JSplitPane;
 
-import bilib.component.PanelImage;
 import bilib.table.CustomizedTable;
 import bilib.tools.NumFormat;
 import deconvolution.Deconvolution;

src/main/java/deconvolution/capsule/PSFCapsule.java

@@ -31,11 +31,8 @@
 
 package deconvolution.capsule;
 
-import java.awt.Dimension;
-
 import javax.swing.JSplitPane;
 
-import bilib.component.PanelImage;
 import bilib.table.CustomizedTable;
 import bilib.tools.NumFormat;
 import deconvolution.Deconvolution;

src/main/java/deconvolution/capsule/ReportCapsule.java

@@ -31,11 +31,8 @@
 
 package deconvolution.capsule;
 
-import java.awt.Dimension;
-
 import javax.swing.JSplitPane;
 
-import bilib.component.PanelImage;
 import bilib.table.CustomizedTable;
 import deconvolution.Deconvolution;
 import signal.RealSignal;

src/main/java/deconvolution/capsule/ResourcesCapsule.java

@@ -35,14 +35,12 @@
 import java.awt.GridLayout;
 import java.awt.event.ActionEvent;
 import java.awt.event.ActionListener;
-import java.awt.event.MouseEvent;
 import java.util.ArrayList;
 
 import javax.swing.JButton;
 import javax.swing.JPanel;
 import javax.swing.JSplitPane;
 
-import bilib.table.CustomizedTable;
 import deconvolution.Deconvolution;
 import deconvolutionlab.system.AbstractMeter;
 import deconvolutionlab.system.FFTMeter;