Commit c8e09488 authored by Andrey Filippov's avatar Andrey Filippov
Browse files

Added bicubic interpolation

parent 42058891

src/main/java/com/elphel/imagej/cameras/CLTParameters.java

+6 −0
Original line number Diff line number Diff line
@@ -493,6 +493,7 @@ public class CLTParameters { 
	public boolean    gmap_render_hdr =      true; //  generate textures w/o normalization to generate undistorted

	public boolean    gmap_en =              true; //  generate ground map from a drone (enables gmap_render_hdr)

	public boolean    gmap_parallel_proj =   true; // Use parallel projection (map)

	public boolean    gmap_bicubic =         true; // Use bicubic interpolation (false - bilinear)

	public boolean    gmap_update_range =    false; // for parallel only

	

	// extracting ground projection plane
@@ -1674,6 +1675,7 @@ public class CLTParameters { 
		properties.setProperty(prefix+"gmap_render_hdr",            this.gmap_render_hdr+"");            // boolean 

		properties.setProperty(prefix+"gmap_en",                    this.gmap_en+"");                    // boolean

		properties.setProperty(prefix+"gmap_parallel_proj",         this.gmap_parallel_proj+"");         // boolean

		properties.setProperty(prefix+"gmap_bicubic",               this.gmap_bicubic+"");               // boolean

		properties.setProperty(prefix+"gmap_update_range",          this.gmap_update_range+"");          // boolean

		

		properties.setProperty(prefix+"gmap_use_lma",               this.gmap_use_lma+"");               // boolean
@@ -2726,6 +2728,7 @@ public class CLTParameters { 
		if (properties.getProperty(prefix+"gmap_render_hdr")!=null)      this.gmap_render_hdr=Boolean.parseBoolean(properties.getProperty(prefix+   "gmap_render_hdr"));

		if (properties.getProperty(prefix+"gmap_en")!=null)              this.gmap_en=Boolean.parseBoolean(properties.getProperty(prefix+           "gmap_en"));

		if (properties.getProperty(prefix+"gmap_parallel_proj")!=null)   this.gmap_parallel_proj=Boolean.parseBoolean(properties.getProperty(prefix+"gmap_parallel_proj"));

		if (properties.getProperty(prefix+"gmap_bicubic")!=null)         this.gmap_bicubic=Boolean.parseBoolean(properties.getProperty(prefix+      "gmap_bicubic"));

		if (properties.getProperty(prefix+"gmap_update_range")!=null)    this.gmap_update_range=Boolean.parseBoolean(properties.getProperty(prefix+ "gmap_update_range"));

		if (properties.getProperty(prefix+"gmap_use_lma")!=null)         this.gmap_use_lma=Boolean.parseBoolean(properties.getProperty(prefix+      "gmap_use_lma"));

		if (properties.getProperty(prefix+"gmap_discard_low")!=null)     this.gmap_discard_low=Double.parseDouble(properties.getProperty(prefix+    "gmap_discard_low"));
@@ -4031,6 +4034,8 @@ public class CLTParameters { 
				"Render fixed-scale projection of the 3D model to a ground plane, such as for the UAS-generated imagery.");

		gd.addCheckbox ("Parallel (not center) projection (maps)", this.gmap_parallel_proj, // true; // enable change FG pixel to opaque from transparent

				"Parallel-project objects to a plane surface. If unchecked - use center (from the camera) projection.");

		gd.addCheckbox ("Use bi-cubic texture interpolation",      this.gmap_bicubic, // true; // enable change FG pixel to opaque from transparent

				"Use bi-cubic texture interpolation. False - old way, is bilinear.");

		gd.addCheckbox ("Fit all model (parallel only)",           this.gmap_update_range, // true; // enable change FG pixel to opaque from transparent

				"Recalculate (increase) image size to fit all model elements. Unchecked - limit by plane intersection with the camera FOV.");
@@ -5307,6 +5312,7 @@ public class CLTParameters { 
		this.gmap_render_hdr=       gd.getNextBoolean();

		this.gmap_en=               gd.getNextBoolean();

		this.gmap_parallel_proj=    gd.getNextBoolean();

		this.gmap_bicubic=          gd.getNextBoolean();

		this.gmap_update_range=     gd.getNextBoolean();

		this.gmap_use_lma=          gd.getNextBoolean();

		this.gmap_discard_low=      gd.getNextNumber();

src/main/java/com/elphel/imagej/tileprocessor/TexturedModel.java

+14 −7
Original line number Diff line number Diff line
@@ -2458,9 +2458,9 @@ public class TexturedModel { 
	{

		final boolean map_en =            clt_parameters.gmap_en;

		final boolean render_hdr =        clt_parameters.gmap_render_hdr || map_en;// true; //false; // true; // generate textures w/o normalization to generate undistorted

		final boolean use_parallel_proj = clt_parameters.gmap_parallel_proj; // Use parallel projection (map)

		final boolean use_parallel_proj = clt_parameters.gmap_parallel_proj; // Use parallel projection (map)\

		final boolean bicubic =           clt_parameters.gmap_bicubic;       // Use bicubic interpolation (false - bilinear)

		final boolean update_range =      clt_parameters.gmap_update_range; // for parallel only

		

		final boolean   use_lma =         clt_parameters.gmap_use_lma ; // true;  // ;

		final double    discard_low =     clt_parameters.gmap_discard_low ; //0.01;  // fraction of all pixels

		final double    discard_high =    clt_parameters.gmap_discard_high ; //0.5;   // fraction of all pixels
@@ -3178,11 +3178,17 @@ public class TexturedModel { 
					boolean last_is_alpha = true; // last channel in textures slices is alpha

					double [][] full_render_z;

					if (use_parallel_proj) {

						if (bicubic) {

							full_render_z =render3D.render3dPlaneParallelProjBiCubic(

									tri_meshes,            // final ArrayList<TriMesh> tri_meshes,

									last_is_alpha,         // final boolean            last_is_alpha,

									debugLevel);           //int                      debugLevel)

						} else {

							full_render_z =render3D.render3dPlaneParallelProj(

									tri_meshes,            // final ArrayList<TriMesh> tri_meshes,

									last_is_alpha,         // final boolean            last_is_alpha,

								//							scenes[ref_index],     //final QuadCLT            ref_scene,    // all coordinates relative to this scene

									debugLevel);           //int                      debugLevel)

						}

					} else {

						full_render_z =render3D.render3dPlaneCenterProj(

								tri_meshes,            // final ArrayList<TriMesh> tri_meshes,
@@ -3193,6 +3199,7 @@ public class TexturedModel { 
					//			String model_name = ref_scene.correctionsParameters.getModelName(ref_scene.getImageName());

					String suffix ="-RECT"; 

					suffix +="-PIX"+pix_size * hdr_whs[2];

					suffix += bicubic ? "-BC":"-BL";

					if (gsmth_enable) {

						suffix +=mixed_flat ? "-FLAT_MIX":"-FLAT_CLN"; // flattened ground - mixed or clean

					}

src/main/java/com/elphel/imagej/x3d/export/Render3D.java

+240 −0
Original line number Diff line number Diff line
@@ -26,6 +26,7 @@ import java.util.ArrayList; 
import java.util.Arrays;

import java.util.concurrent.atomic.AtomicInteger;



import org.apache.commons.math3.analysis.interpolation.PiecewiseBicubicSplineInterpolatingFunction;

import org.apache.commons.math3.geometry.euclidean.threed.Vector3D;

import org.apache.commons.math3.util.FastMath;
@@ -446,6 +447,245 @@ public class Render3D { 
		ImageDtt.startAndJoin(threads);

		return full_rendered;

	}

//https://commons.apache.org/proper/commons-math/javadocs/api-3.6.1/org/apache/commons/math3/analysis/interpolation/PiecewiseBicubicSplineInterpolatingFunction.html	

	public  double [][] render3dPlaneParallelProjBiCubic(

			final ArrayList<TriMesh> tri_meshes,

			final boolean            last_is_alpha,

			int                      debugLevel){  // debug level

		// TODO: add crop - add to the caller

		if ((tri_meshes == null) || tri_meshes.isEmpty() || (tri_meshes.get(0).getTexturePixels() == null)) {

			return null;

		}

		final boolean export_z = true;

		final int dbg_ipix=1673752;

		// get total number of triangles

		int num_tri=0;

		for (TriMesh tri: tri_meshes) {

			num_tri += tri.getTriangles().length;

		}

		int indx=0;

		int num_mesh = 0;

		final int z_index = tri_meshes.get(0).getTexturePixels().length;

		final int [][] tri_index = new int[num_tri][2];

		for (TriMesh tri: tri_meshes) {

			int num_tri_mesh = tri.getTriangles().length;

			for (int i = 0; i < num_tri_mesh; i++)

			{

				tri_index[indx][0] =   num_mesh;

				tri_index[indx++][1] = i;

			}

			num_mesh++;

		}

		final PiecewiseBicubicSplineInterpolatingFunction [][] pbsif =

				new PiecewiseBicubicSplineInterpolatingFunction[tri_meshes.size()][z_index];

		final double [][] x = new double [tri_meshes.size()][];

		final double [][] y = new double [tri_meshes.size()][];

		final double [][][][] f = new double [tri_meshes.size()][z_index][][];

		

		for (int nm =0; nm < pbsif.length; nm++) {

			TriMesh tri = tri_meshes.get(nm);

			int texture_width=  tri.getTextureWidth();

			int texture_height= tri.getTextureHeight();

			double [][] texture =       tri.getTexturePixels();

			x[nm] = new double[texture_width];

			y[nm] = new double[texture_height];

			for (int i = 0; i < x[nm].length; i++) {

				x[nm][i] = i;  // add 0.5? Linear does not

			}

			for (int i = 0; i < y[nm].length; i++) {

				y[nm][i] = i;  // add 0.5? Linear does not

			}

			for (int chn = 0; chn < z_index; chn++) {

				f[nm][chn] = new double [texture_width][texture_height];

				for (int row = 0; row < texture_height; row++) {

					for (int col = 0; col < texture_width; col++) {

						f[nm][chn][col][row] = texture[chn][col + row * texture_width];	

					}

					/*

					System.arraycopy(

							texture[chn],

							row * texture_width,

							f[nm][chn][row],

							0,

							texture_width);

							*/

				}

				// clones x,y,f (but not f[]). f is in columns, not rows! f[width][height]

				pbsif[nm][chn] = new PiecewiseBicubicSplineInterpolatingFunction(x[nm],y[nm],f[nm][chn]);

			}

		}		

		

		

		if (debugLevel > -2) {

			System.out.println("Prepare to render "+num_tri+" triangles in "+num_mesh+" meshes");

		}

		final double [][] full_rendered = new double[z_index+ (export_z? 1:0)][out_width * out_height];

		int alpha_index = last_is_alpha ? (z_index - 1) : z_index;

		for (int chn = 0; chn < alpha_index; chn++) {

			Arrays.fill(full_rendered[chn], Double.NaN);

		}

		// create z-buffer array per each thread, in the end - merge them

		final Thread[] threads = ImageDtt.newThreadArray(THREADS_MAX);

		final AtomicInteger ai = new AtomicInteger(0);

		final AtomicInteger ati = new AtomicInteger(0);

		final double [][][] rendered = new double [threads.length][full_rendered[0].length][];

		for (int ithread = 0; ithread < threads.length; ithread++) {

			threads[ithread] = new Thread() {

				public void run() {

					int ti = ati.getAndIncrement();

					double [][] rend = rendered[ti]; // this thread rendered results

					for (int indx = ai.getAndIncrement(); indx < tri_index.length; indx = ai.getAndIncrement()) {

						TriMesh tri=tri_meshes.get(tri_index[indx][0]); // mesh to process;

						int tri_indx = tri_index[indx][1]; // triangle index in a mesh

						double [][] texture =       tri.getTexturePixels();

						int texture_width=          tri.getTextureWidth();

						int texture_height=         tri.getTextureHeight();

						double[][] mesh_coord =     tri.getCoordinates();

						double[][] mesh_tex_coord = tri.getTexCoord();

						int [] triangle = tri.getTriangles()[tri_indx];

						double [][] tri_out2=  new double[3][];

						double [][] tri_text2= new double[3][];

						double [][] min_max_xyz = new double[3][2];

						for (int i = 0; i < 3; i++) {

							double [] gxyz=ErsCorrection.applyXYZATR(toground, mesh_coord[triangle[i]]);

							tri_out2[i] =  new double [] {

									pixel_per_m * (gxyz[0] - x0_y0[0]), 

									pixel_per_m * (gxyz[1] - x0_y0[1]),

									-pixel_per_m * gxyz[2] // make it positive?  

							};

							//[2] - distance from the camera in "pixels" - same linear scale as on the ground. lower values obscure higher.

							tri_text2[i] = mesh_tex_coord[triangle[i]];

							for (int j = 0; j < 3; j++) {

								if ((i==0) || (tri_out2[i][j] < min_max_xyz[j][0])) min_max_xyz[j][0] = tri_out2[i][j]; 

								if ((i==0) || (tri_out2[i][j] > min_max_xyz[j][1])) min_max_xyz[j][1] = tri_out2[i][j]; 

							}

						}

						// Check plane direction

						double [] d01 = new double [] {tri_out2[1][0]-tri_out2[0][0], tri_out2[1][1]-tri_out2[0][1]};

						double [] d02 = new double [] {tri_out2[2][0]-tri_out2[0][0], tri_out2[2][1]-tri_out2[0][1]};

						if (!cross2ccw(d02,d01)) {

							continue;

						}

//						PiecewiseBicubicSplineInterpolatingFunction pbs = pbsf[tri_index[indx][0]]

						int ipx_min = (int) Math.floor(min_max_xyz[0][0]);

						int ipx_max = (int) Math.ceil (min_max_xyz[0][1]);

						int ipy_min = (int) Math.floor(min_max_xyz[1][0]);

						int ipy_max = (int) Math.ceil (min_max_xyz[1][1]);

						// apply bounds

						if (ipx_min < 0) ipx_min = 0;

						if (ipy_min < 0) ipy_min = 0;

						if (ipx_max >= out_width) ipx_max = out_width - 1;

						if (ipy_max >= out_height) ipy_max = out_height - 1;

						if ((ipx_min > ipx_max) || (ipy_min > ipy_max)) {

							continue; // triangle completely outside rendering are

						}

						

						//pbsif[nm][chn]

						// vector from 1 to 2

						double [] t01 = new double [] {tri_text2[1][0]-tri_text2[0][0], tri_text2[1][1]-tri_text2[0][1]};

						double [] t02 = new double [] {tri_text2[2][0]-tri_text2[0][0], tri_text2[2][1]-tri_text2[0][1]};

						double [][] orto2 = orthonorm2(d01, d02);

						double [] d12 = new double [] {tri_out2[2][0]-tri_out2[1][0], tri_out2[2][1]-tri_out2[1][1]};

						for (int ipy = ipy_min; ipy <= ipy_max; ipy++) {

							for (int ipx = ipx_min; ipx <= ipx_max; ipx++) {

								// check it is inside triangle

								double [] d0p = new double[] {ipx-tri_out2[0][0],ipy-tri_out2[0][1]};

								if (!cross2ccw(d0p,d01)) continue; 

								if (!cross2ccw(d02,d0p)) continue; 

								double [] d1p = new double[] {ipx-tri_out2[1][0],ipy-tri_out2[1][1]};

								if (!cross2ccw(d1p,d12)) continue;

								int ipix = ipx +(out_height - 1 -ipy) * out_width; // Y goes down

								if (ipix== dbg_ipix) {

									System.out.println("ipix="+ipix);

								}

								// See if the rendered pixel is closer than the closest of the corners

								if ((rend[ipix] != null ) && (rend[ipix][z_index] < min_max_xyz[2][0])) {

									continue; 

								}

								double kx = dot2(d0p, orto2[0]);

								double ky = dot2(d0p, orto2[1]);

								// interpolate z

								double z_interp = tri_out2[0][2] +

										kx * (tri_out2[1][2]-tri_out2[0][2]) +

										ky * (tri_out2[2][2]-tri_out2[0][2]);

								if ((rend[ipix] != null ) && (rend[ipix][z_index] < z_interp)) {

									continue;

								}

								// Get corresponding texture coordinates

								double text_x = tri_text2[0][0] + kx * t01[0] + ky*t02[0]; // texture relative coordinates (0,1)

								double text_y = tri_text2[0][1] + kx * t01[1] + ky*t02[1]; // y - up!

								double px = text_x * texture_width - 0.5; // (0.0,0.0) - center of top-left texture pixel

								double py = (1.0-text_y) * texture_height -0.5;

								int ipx0 = (int) Math.floor(px);

								int ipy0 = (int) Math.floor(py);

								double fx = px - ipx0;

								double fy = py - ipy0;

								int ipx1 = ipx0+1;

								int ipy1 = ipy0+1;

								if ((ipx1 < 0) || (ipy1 < 0) || (ipx0 >= texture_width) || (ipy0 >= texture_width)) {

									continue; // outside bounds

								}

								// limit if just on the edge

								if (ipx0 < 0) ipx0=ipx1;

								if (ipy0 < 0) ipy0=ipy1;

								if (ipx1 >= texture_width)  ipx1=ipx0;

								if (ipy1 >= texture_height) ipy1=ipy0;

								int indx00 = ipx0+texture_width*ipy0;

								int indx10 = ipx1+texture_width*ipy0;

								int indx01 = ipx0+texture_width*ipy1;

								int indx11 = ipx1+texture_width*ipy1;

								double [] pix_val = new double[z_index+1];

								pix_val[z_index] = z_interp;

								for (int chn = 0; chn < z_index; chn++) { // just for testing - use bilinear first

									pix_val[chn] =

											(1.0 - fy) * (1.0 - fx) * texture[chn][indx00] + 

											(1.0 - fy) * (      fx) * texture[chn][indx10] + 

											(      fy) * (1.0 - fx) * texture[chn][indx01] + 

											(      fy) * (      fx) * texture[chn][indx11];

									if (pbsif[tri_index[indx][0]][chn].isValidPoint(px, py)) { // tghen overwrite with bicubic

										pix_val[chn] = pbsif[tri_index[indx][0]][chn].value(px,py);   

									}

								}

								// handle alpha

								if (last_is_alpha && (pix_val[z_index-1] < 0.5)) {

									continue; // low alpha -> transparent

								}

								rend[ipix] = pix_val;

							}

						}

						// min_max_xyz[2]

						//num_col_chn

						// projectToPlanePixels

						// getCoordinates()

						//getTexCoord()

					}

				}

			};

		}		      

		ImageDtt.startAndJoin(threads);

		// merge partial renders:

		ai.set(0);

		for (int ithread = 0; ithread < threads.length; ithread++) {

			threads[ithread] = new Thread() {

				public void run() {

					for (int indx = ai.getAndIncrement(); indx < full_rendered[0].length; indx = ai.getAndIncrement()) {

						double z = Double.NaN;

						for (int sub_render = 0; sub_render < rendered.length; sub_render++) if (rendered[sub_render][indx] != null){

							if (!(rendered[sub_render][indx][z_index] <= z)) { // OK previous NaN

								z = rendered[sub_render][indx][z_index];

								for (int chn = 0; chn < full_rendered.length; chn++) { // z_index; chn++) {

									full_rendered[chn][indx] = rendered[sub_render][indx][chn];

								}

							}

						}

					}

				}

			};

		}		      

		ImageDtt.startAndJoin(threads);

		return full_rendered;

	}

	

	

	public  double [][] render3dPlaneCenterProj(

			final ArrayList<TriMesh> tri_meshes,