package com.elphel.imagej.calibration;
/**
** -----------------------------------------------------------------------------**
** PixelMapping.java
**
** Using camera calibration files to resample images for equirectangular projection
**
**
** Copyright (C) 2012 Elphel, Inc.
**
** -----------------------------------------------------------------------------**
**
** PixelMapping.java is free software: you can redistribute it and/or modify
** it under the terms of the GNU General Public License as published by
** the Free Software Foundation, either version 3 of the License, or
** (at your option) any later version.
**
** This program is distributed in the hope that it will be useful,
** but WITHOUT ANY WARRANTY; without even the implied warranty of
** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
** GNU General Public License for more details.
**
** You should have received a copy of the GNU General Public License
** along with this program. If not, see .
** -----------------------------------------------------------------------------**
**
*/
import java.awt.Point;
import java.awt.Rectangle;
import java.io.IOException;
import java.util.ArrayList;
import java.util.Arrays;
import java.util.BitSet;
import java.util.Comparator;
import java.util.List;
import java.util.concurrent.atomic.AtomicInteger;
import javax.swing.SwingUtilities;
import com.elphel.imagej.cameras.EyesisCorrectionParameters;
import com.elphel.imagej.common.DoubleFHT;
import com.elphel.imagej.common.DoubleGaussianBlur;
import com.elphel.imagej.common.ShowDoubleFloatArrays;
import com.elphel.imagej.jp4.JP46_Reader_camera;
import com.elphel.imagej.readers.EyesisTiff;
import Jama.Matrix;
import ij.IJ;
import ij.ImagePlus;
import ij.ImageStack;
import ij.Prefs;
//import ij.gui.GenericDialog;
import ij.io.FileSaver;
import ij.io.Opener;
import ij.process.ColorProcessor;
import ij.text.TextWindow;
import loci.common.services.DependencyException;
import loci.common.services.ServiceException;
import loci.formats.FormatException;
public class PixelMapping {
public SensorData [] sensors;
public int panoWidth=0;
public int panoHeight=0;
public double panoDegreesPerPixel=Double.NaN; // needs to be initialized
public double panoLongitudeLeft=Double.NaN;
public double panoLongitudeRight=Double.NaN;
public double panoLatitudeTop=Double.NaN;
public double panoLatitudeBottom=Double.NaN;
public int debugLevel=1;
public int lanczosA=3;
public int oversampled=2;
public int binsPerHalfPixel=50;
public double [][][][] lanczos=null;
public int maxSensors=100;
InterSensor lastUsedInterSensor=null;
// public enum cellTypes {EMPTY_CELL, BAD_CELL, OLD_CELL,NEW_CELL}
public PixelMapping (String defaultPath, int debugLevel){
this.debugLevel=debugLevel;
String [] extensions={".calib-tiff"};
String [] defaultPaths=((defaultPath==null) || defaultPath.equals(""))?null:(new String[1]);
if (defaultPaths!=null) defaultPaths[0]=defaultPath;
MultipleExtensionsFileFilter parFilter = new MultipleExtensionsFileFilter("",extensions,"distortion calibration .calib-tiff files");
String [] calibFiles=CalibrationFileManagement.selectFiles(false,
"Select camera sub-modules calibration files",
"Select",
parFilter,
defaultPaths); // String [] defaultPaths);
if ((calibFiles==null) || (calibFiles.length==0)) {
IJ.showMessage("No files selected");
}
int maxChannel=0;
for (int i=0;imaxChannel) maxChannel=channel;
}
this.sensors=new SensorData[maxChannel+1];
for (int i=0;i 0)? num_channels: this.maxSensors];
for (int i=0;i=16) && (channel < 16)) {
if (channel_file != channel) {
System.out.println ("****Temporary fix to use old sensor files, COMMENT OUT when done debugging****");
System.out.println ("i="+i+", channel_file="+channel_file+", channel ="+channel+" -> "+channel_file);
channel = channel_file; // i;
}
if ((channel < first_channel) || ((num_channels > 0) && (channel >= (first_channel + num_channels)))) {
continue; // wrong channels
}
channel -= first_channel;
if (update_channel) {
sensorData.setChannel(channel);
}
this.sensors[channel]=sensorData;
if (channel>maxChannel) maxChannel=channel;
}
SensorData [] tmpSensors=this.sensors.clone();
this.sensors=new SensorData[maxChannel+1];
for (int i=0;i=this.sensors.length)) return null;
return this.sensors[channel].getPath();
}
public int getSubChannel(int channel ){
if ((channel<0) || (channel>=this.sensors.length)) {
System.out.println("ERROR: getSubChannel("+channel+"): channel out of range");
return -1;
}
// System.out.println("getSubChannel("+channel+") => "+this.sensors[channel].getSubChannel());
return this.sensors[channel].getSubChannel();
}
public int getSubChannelSilent(int channel ){
if ((channel<0) || (channel>=this.sensors.length)) {
return -1;
}
// System.out.println("getSubChannel("+channel+") => "+this.sensors[channel].getSubChannel());
return this.sensors[channel].getSubChannel();
}
public int getSubCamera(int channel ){
if ((channel<0) || (channel>=this.sensors.length)) return -1;
return this.sensors[channel].getSubCamera();
}
// Enumerated calibration file
public boolean isChannelAvailable(int channel){
return (this.sensors != null) && (channel>=0) && (channel= 0)) {
return true;
}
return false;
}
// if subcameras are not used (-1), as in Quads, return {channel index} or null if alien file
public int [] getSensorWH(int chn) {
if ((sensors == null) || (chn < 0) || (chn >= sensors.length)) {
return null;
}
return sensors[chn].getSensorWH();
}
public int [] channelsForSubCamera(int subCamera){
if (subCamera < 0) {
return null;
}
if (this.sensors == null) return null;
if (this.debugLevel>1) {
System.out.print("channelsForSubCamera("+subCamera+"),this.sensors.length="+this.sensors.length+": ");
}
if (!subcamerasUsed()) {
for (int i=0;i cam_port = new ArrayList();
for (int i=0;i() {
@Override
public int compare(Point o1, Point o2) {
return (o1.x>o2.x)? 1:((o1.x < o2.x)?-1:(o1.y > o2.y)? 1:((o1.y < o2.y)?-1:0));
}
});
// debugging:
if (subCamera >= cam_port_arr.length) {
if (this.debugLevel>1) {
System.out.println("Error: Subcamera "+subCamera+" > that total number of sensor ports in the system = "+cam_port_arr.length);
}
return null;
}
if (this.debugLevel>1) {
System.out.println("----- This filename subcamera "+subCamera+": physical camera "+cam_port_arr[subCamera].x+", sensor_port "+cam_port_arr[subCamera].y);
} else if (this.debugLevel > -2) {
System.out.print(".");
}
int numChannels=0;
for (int i=0;i=0) && (channel=this.sensors.length))return null;
return this.sensors[channel].getBayerFlatFieldFloat(
bayer,
range);
}
public double [] getBayerFlatField(
int channel,
int [][] bayer){ //{{1,0},{2,1}} GR/BG
if ((this.sensors == null) || (channel<0) && (channel>=this.sensors.length))return null;
return this.sensors[channel].getBayerFlatField(
// width,
// height,
bayer);
}
public int [][] getDefectsXY(
int channel){
if ((this.sensors == null) || (channel<0) && (channel>=this.sensors.length))return null;
return this.sensors[channel].getDefectsXY();
}
public double[] getDefectsDiff(
int channel){
if ((this.sensors == null) || (channel<0) && (channel>=this.sensors.length))return null;
return this.sensors[channel].getDefectsDiff();
}
///SensorData
/*
public float [] getBayerFlatFieldFloat(
int width,
int height,
int [][] bayer){ //{{1,0},{2,1}} GR/BG
*
int maxChannel=0;
for (int i=0;imaxChannel) maxChannel=channel;
}
this.sensors=new SensorData[maxChannel+1];
for (int i=0;imaxChannel) maxChannel=channel;
}
if (this.sensors==null) {
this.sensors=new SensorData[maxChannel+1];
for (int i=0;i0) System.out.println(msg);
if (this.sensors[channel]==null) this.sensors[channel] =new SensorData(channelPath,true);
else this.sensors[channel].createEquirectangularMap(channelPath);
}
}
public void loadChannelEquirectangularMap(
int channel,
String path){
if ((this.sensors==null) || (this.sensors[channel]==null)){
String msg="Sensor "+channel+" data is not initialized";
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
// this.sensors[channel] =new SensorData(path,true);
this.sensors[channel].createEquirectangularMap(path);
}
public void resampleToEquirectangular(
String [] paths,
int sourceImageScale,
boolean showResults,
boolean saveResults,
int maxThreads){
for (int i=0;ichannel) && (this.sensors[channel]!=null)){
this.sensors[channel].equirectangularMap=null;
}
}
public void deleteEquirectangularMapFull(int channel){
if ((this.sensors!=null) && (this.sensors.length>channel) && (this.sensors[channel]!=null) && (this.sensors[channel].equirectangularMap!=null)){
this.sensors[channel].equirectangularMap.map=null;
}
}
public boolean isEquirectangularMapAvailable(int channel){
if ((this.sensors!=null) && (this.sensors.length>channel) && (this.sensors[channel]!=null)){
return((this.sensors[channel].equirectangularMap!=null) && (this.sensors[channel].equirectangularMap.partialMap!=null));
}
return false;
}
public boolean isPlaneMapMapAvailable(int channel){
if ((this.sensors!=null) && (this.sensors.length>channel) && (this.sensors[channel]!=null)){
return((this.sensors[channel].interSensor!=null));
}
return false;
}
public ImagePlus resampleToEquirectangular(
ImagePlus imp,
int channel,
int sourceImageScale,
int maxThreads){
if (imp==null) return null;
if ((this.sensors==null) || (this.sensors.length<(channel+1) || (this.sensors[channel]==null))){
IJ.showMessage("No sensor data for channel "+channel);
return null;
}
if (( this.sensors[channel].equirectangularMap==null) || ( this.sensors[channel].equirectangularMap.partialMap==null)){
IJ.showMessage("No equirectangular map for channel "+channel);
return null;
}
if (this.lanczos==null) generateLanczosStack();
if (imp.getType()==ImagePlus.COLOR_RGB) {
return resampleToEquirectangularRGB24(imp, channel, sourceImageScale, maxThreads);
} else if (imp.getStackSize()>=3) {
return resampleToEquirectangularRGBFP32(imp, channel, sourceImageScale, maxThreads);
}
IJ.showMessage("Not yet implemented for this image type");
return null;
}
/**
* Warp 3-slice 32-bit FP stack (RGB) to 4-slice stack RGBA
* @param imp 3-slice image stack with straight image
* @param channel sensor channel number
* @param scale source image scale (normally 2x)
* @param maxThreads
* @return 4-slice RGBA image, floating point 32 bits per sample per color channel
*/
public ImagePlus resampleToEquirectangularRGBFP32(
ImagePlus imp,
int channel,
final int scale,
int maxThreads){
final SensorData.EquirectangularMap erm=this.sensors[channel].equirectangularMap;
final int width=imp.getWidth();
final int height=imp.getHeight();
if (scale==0) {
IJ.showMessage("resampleToEquirectangularRGBFP32(): Image has less resolution than the map, can not proceed");
return null;
}
if (imp.getStackSize()<3) {
IJ.showMessage("resampleToEquirectangularRGBFP32(): Expecting 3-layer stack(R,G,B)");
return null;
}
final float [][] imagePixels= new float [3][];
for (int c=0;c0){ // do not convolve pixels with alpha=0;
double x=scale*erm.partialMap[0][oIndex];
double y=scale*erm.partialMap[1][oIndex];
int ix= (int) Math.round(x);
int iy= (int) Math.round(y);
double dx=x-ix;
double dy=y-iy;
int indxX= (int) Math.round((2*dx+1)*binsPerHalfPixel);
int indxY= (int) Math.round((2*dy+1)*binsPerHalfPixel);
double [][] lk=lanczos[indxY][indxX];
for (int i=0;i=height) ipy=height-1;
int baseI=ipy*width;
for (int j=0;j=width) ipx=width-1;
/*
int pix=imagePixels[baseI+ipx];
RGBA[1]+=((pix>>16) & 0xff)*lk[i][j]; // R
RGBA[2]+=((pix>>8) & 0xff)*lk[i][j]; // G
RGBA[3]+=( pix & 0xff)*lk[i][j]; // B
*/
if ((baseI+ipx)>imagePixels[0].length){
System.out.println("baseI+ipx="+(baseI+ipx)+" baseI="+baseI+" ipx="+ipx+" i="+i+" j="+j+
" erm.mapWOI.width="+erm.mapWOI.width+" erm.mapWOI.height="+erm.mapWOI.height);
}
RGBA[1]+=imagePixels[0][baseI+ipx]*lk[i][j]; // R
RGBA[2]+=imagePixels[1][baseI+ipx]*lk[i][j]; // G
RGBA[3]+=imagePixels[2][baseI+ipx]*lk[i][j]; // B
}
}
outPixels[0][oIndex]=(float) RGBA[1];
outPixels[1][oIndex]=(float) RGBA[2];
outPixels[2][oIndex]=(float) RGBA[3];
/*
int c=0;
for (int i=0;i<4;i++){
if (RGBA[i]< 0.0) RGBA[i]=0.0;
else if (RGBA[i] > 255.0) RGBA[i]=255.0;
c=(c<<8) | ((int) RGBA[i]);
outPixels [oIndex]=c;
}
*/
} else {
outPixels[0][oIndex]=0.0F;
outPixels[1][oIndex]=0.0F;
outPixels[2][oIndex]=0.0F;
// outPixels[oIndex]=0;
}
}
final int numFinished=opyDoneAtomic.getAndIncrement();
// IJ.showProgress(progressValues[numFinished]);
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
// Here, we can safely update the GUI
// because we'll be called from the
// event dispatch thread
IJ.showProgress(numFinished,erm.mapWOI.height);
}
});
}
}
};
}
startAndJoin(threads);
ImageStack outStack=new ImageStack(erm.mapWOI.width,erm.mapWOI.height);
outStack.addSlice("Red", outPixels[0]);
outStack.addSlice("Green", outPixels[1]);
outStack.addSlice("Blue", outPixels[2]);
outStack.addSlice("Alpha", outPixels[3]);
// ColorProcessor cp=new ColorProcessor(erm.mapWOI.width,erm.mapWOI.height);
// cp.setPixels(outPixels);
ImagePlus impOut=new ImagePlus(imp.getTitle()+"_EQR",outStack);
impOut.setProperty("channel", ""+erm.channel); // image channel
impOut.setProperty("XPosition", ""+erm.mapWOI.x); // real XPosition as tiff tag 0x11e is rational (5)
impOut.setProperty("YPosition", ""+erm.mapWOI.y); // real XPosition as tiff tag 0x11f is rational (5)
impOut.setProperty("ImageFullWidth", "" +erm.pixelsHorizontal);// real ImageFullWidth as tiff tag 0x8214 is rational (5)
impOut.setProperty("ImageFullLength", ""+erm.pixelsVertical); // real ImageFullLength as tiff tag 0x8215 is rational (5)
(new JP46_Reader_camera(false)).encodeProperiesToInfo(impOut);
impOut.getProcessor().resetMinAndMax();
return impOut;
}
public ImagePlus resampleToEquirectangularRGB24(
ImagePlus imp,
int channel,
final int scale,
int maxThreads){
final SensorData.EquirectangularMap erm=this.sensors[channel].equirectangularMap;
final int width=imp.getWidth();
final int height=imp.getHeight();
if (scale==0) {
IJ.showMessage("Image has less resolution than the map, can not proceed");
return null;
}
final int [] imagePixels=(int []) imp.getProcessor().getPixels();
final int [] outPixels=new int [erm.mapWOI.width*erm.mapWOI.height];
if (this.oversampled!=scale) generateLanczosStack(
this.lanczosA,
scale,
this.binsPerHalfPixel);
final double [][][][] lanczos=this.lanczos;
final int center=this.lanczos[0][0][0].length/2;
final Thread[] threads = newThreadArray(maxThreads);
final AtomicInteger opyAtomic = new AtomicInteger(0);
final AtomicInteger opyDoneAtomic = new AtomicInteger(1);
final int binsPerHalfPixel=this.binsPerHalfPixel;
// for (int opy=0;opy0){ // do not convolve pixels with alpha=0;
double x=scale*erm.partialMap[0][oIndex];
double y=scale*erm.partialMap[1][oIndex];
int ix= (int) Math.round(x);
int iy= (int) Math.round(y);
double dx=x-ix;
double dy=y-iy;
int indxX= (int) Math.round((2*dx+1)*binsPerHalfPixel);
int indxY= (int) Math.round((2*dy+1)*binsPerHalfPixel);
double [][] lk=lanczos[indxY][indxX];
for (int i=0;i=height) ipy=height-1;
int baseI=ipy*width;
for (int j=0;j=width) ipx=width-1;
int pix=imagePixels[baseI+ipx];
RGBA[1]+=((pix>>16) & 0xff)*lk[i][j]; // R
RGBA[2]+=((pix>>8) & 0xff)*lk[i][j]; // G
RGBA[3]+=( pix & 0xff)*lk[i][j]; // B
}
}
int c=0;
for (int i=0;i<4;i++){
if (RGBA[i]< 0.0) RGBA[i]=0.0;
else if (RGBA[i] > 255.0) RGBA[i]=255.0;
c=(c<<8) | ((int) RGBA[i]);
outPixels [oIndex]=c;
}
} else {
outPixels[oIndex]=0;
}
}
final int numFinished=opyDoneAtomic.getAndIncrement();
// IJ.showProgress(progressValues[numFinished]);
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
// Here, we can safely update the GUI
// because we'll be called from the
// event dispatch thread
IJ.showProgress(numFinished,erm.mapWOI.height);
}
});
}
}
};
}
startAndJoin(threads);
ColorProcessor cp=new ColorProcessor(erm.mapWOI.width,erm.mapWOI.height);
cp.setPixels(outPixels);
ImagePlus impOut=new ImagePlus(imp.getTitle()+"_EQR",cp);
impOut.setProperty("channel", ""+erm.channel); // image channel
impOut.setProperty("XPosition", ""+erm.mapWOI.x); // real XPosition as tiff tag 0x11e is rational (5)
impOut.setProperty("YPosition", ""+erm.mapWOI.y); // real XPosition as tiff tag 0x11f is rational (5)
impOut.setProperty("ImageFullWidth", "" +erm.pixelsHorizontal);// real ImageFullWidth as tiff tag 0x8214 is rational (5)
impOut.setProperty("ImageFullLength", ""+erm.pixelsVertical); // real ImageFullLength as tiff tag 0x8215 is rational (5)
(new JP46_Reader_camera(false)).encodeProperiesToInfo(impOut);
impOut.getProcessor().resetMinAndMax();
return impOut;
}
/**
* Generate a stack of Lanczos kernels to for re-sampling
* @param lacrososA - value of "a" in sinc(x)*sinc(x/a)
* @param oversampled reduce resolution of the image if it was over-sampled (2 for current Eeyesis software)
* @param binsPerHalfPixel calculate for this number of fractional pixels (oversampled pixels, not sensor ones)
* @return [fractional pixels vertical] [fractional pixels horizontal][delta pixel horizontal][ delta pixel vertical]
*/
public double [][][][] generateLanczosStack(){
return generateLanczosStack(
this.lanczosA,
this.oversampled,
this.binsPerHalfPixel);
}
public double [][][][] generateLanczosStack(
int lanczosA,
int oversampled,
int binsPerHalfPixel){
this.lanczosA=lanczosA;
this.oversampled=oversampled;
this.binsPerHalfPixel=binsPerHalfPixel;
if (this.debugLevel>0) System.out.println("Generating Lanczos kernel stack with A="+ this.lanczosA+", oversampled="+oversampled+", binsPerHalfPixel="+binsPerHalfPixel);
double [][][][] lanczos= new double [2*binsPerHalfPixel+1][2*binsPerHalfPixel+1][2*lanczosA*oversampled+1][2*lanczosA*oversampled+1];
double step=0.5/binsPerHalfPixel;
double centerPix=lanczosA*oversampled;
double pi2=Math.PI*Math.PI;
double s=lanczosA/pi2/oversampled;
double sync0=1.0/oversampled;
for (int iy=0;iy=lanczosA)?0.0:((py==0.0)?sync0:(s*Math.sin(py*Math.PI)*Math.sin(py*Math.PI/lanczosA)/(py*py)));
for (int ipx=0;ipx=lanczosA)?0.0:((px==0.0)?sync0:(s*Math.sin(px*Math.PI)*Math.sin(px*Math.PI/lanczosA)/(px*px))));
}
}
}
}
}
this.lanczos=lanczos;
if (this.debugLevel>0) System.out.println("Generated Lanczos kernel stack with A="+ this.lanczosA+", oversampled="+oversampled+", binsPerHalfPixel="+binsPerHalfPixel);
return lanczos;
}
public double [] testLanczosCenter(double [][][][] lanczos){
int size= lanczos.length;
int center= lanczos[0][0].length/2;
double [] testArr=new double [size*size];
for (int i=0;i0) System.out.println("Accumulating channel "+channelNumber);
SensorData.EquirectangularMap erm=this.sensors[channelNumber].equirectangularMap;
if (pixels==null){
this.panoWidth=erm.pixelsHorizontal;
this.panoHeight=erm.pixelsVertical;
pixels= new float [this.panoWidth*this.panoHeight];
for (int i=0;i0){ // alpha
pixels[(iTileLat+woi.y)*this.panoWidth+((iTileLong+woi.x) % this.panoWidth)]+=erm.partialMap[2][tileIndex-1];
}
}
}
}
return pixels;
}
/**
* For up to 32 sensors, calculates in which of them current lat/long is available
* @param threshold - count only pixels with alpha > this threshold (0.0 - all non-zero)
* @return bit masks in scan-line (long, then lat) order
* sets this.panoDegreesPerPixel
*/
public int [] generateOverlapBits( double threshold, int [] sensorList ){
if (sensorList==null){ // all sensors
sensorList=new int [this.sensors.length];
for (int i=0;i0) System.out.println("Accumulating channel "+channelNumber);
SensorData.EquirectangularMap erm=this.sensors[channelNumber].equirectangularMap;
if (usedSensors==null){
this.panoWidth=erm.pixelsHorizontal;
this.panoHeight=erm.pixelsVertical;
usedSensors= new int [this.panoWidth*this.panoHeight];
for (int i=0;ithreshold){ // alpha
usedSensors[(iTileLat+woi.y)*this.panoWidth+((iTileLong+woi.x) % this.panoWidth)] |= (1<=overlapThreshold) numPairs++;
}
int [][] pairs = new int [numPairs][2];
numPairs=0;
for (int i=0;i<(overlap.length-1);i++) for (int j=i+1;j=overlapThreshold) {
pairs[numPairs ][0]=i;
pairs[numPairs++][1]=j;
}
}
return pairs;
}
/*
double threshold=0.01*gd.getNextNumber();
double [][] overlap=PIXEL_MAPPING.overlapPairsAreas(threshold);
PIXEL_MAPPING.listOverlap(overlap);
for (int i=0;i<(overlap.length-1);i++) for (int j=i+1;j=0) indexSuffix--; // remove 1 or 2 digits before period
for (int channelNumber=0;channelNumber0) System.out.println(msg);
sensors[channelNumber].combineDistortionsEquirectangularToSensor(
channelNumber,
longitudeLeft,
longitudeRight,
latitudeTop,
latitudeBottom,
pixelsHorizontal,
imageWidth, //int width,
imageHeight, //int height,
x0, //double x0,
y0, //double y0,
pixelStep, //double pixelStep,
maxThreads);
sensors[channelNumber].equirectangularMap.createPartialMap(longitudeCrop);
if (deleteFull) sensors[channelNumber].equirectangularMap.map=null;
String channelPath=path.substring(0,indexSuffix)+String.format("%02d",channelNumber)+path.substring(indexPeriod);
msg="Saving equirectangular map to "+channelPath;
IJ.showStatus(msg);
if (this.debugLevel>0) System.out.println(msg);
sensors[channelNumber].equirectangularMap.saveMapAsImageStack("EQRCT_MAP"+channelNumber, channelPath);
if (deleteAll) sensors[channelNumber].equirectangularMap=null;
if (this.debugLevel >1) System.out.println("Saved equirectangular map "+channelNumber+" at "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
}
if (this.debugLevel >0) System.out.println("Finished equirectangular maps at "+ IJ.d2s(0.000000001*(System.nanoTime()-startTime),3));
IJ.showStatus("Finished equirectangular maps");
}
public int getNumSubCameras(){
return (this.sensors==null)?0:this.sensors.length;
}
public double [] getParametersVector(int chn){
return this.sensors[chn].getParametersVector();
}
public String getParameterName(int index){
return (new SensorData()).getParameterName(index);
}
public String getParameterDescription(int index){
return (new SensorData()).getParameterDescription(index);
}
public String getParameterUnits(int index){
return (new SensorData()).getParameterUnits(index);
}
public double[][] applyOverlapMap(
String path, //
String imgPathFormat,
String resultDirectory,
int sourceImageScale, // 2.0
boolean saveTiff,
boolean convertToDouble,
int maxThreads,
int debugLevel
){
Opener opener=new Opener();
ImagePlus impMap=opener.openImage("", path);
if (impMap==null) {
String msg="Failed to read inter-sensor overlap map file "+path;
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
if (debugLevel>0) System.out.println("Read "+path+" as an inter-sensor overlap map");
(new JP46_Reader_camera(false)).decodeProperiesFromInfo(impMap);
InterSensor interSensor=new InterSensor(impMap,debugLevel);
this.lastUsedInterSensor=interSensor;
double [][] aYCbCr=convertToDouble?(new double [8][]):null;
for (int iChn=0;iChn<2;iChn++){
int channel=interSensor.channel[iChn];
String sourcePath=String.format(imgPathFormat,channel);
ImagePlus impSrc=opener.openImage("", sourcePath);
if (impSrc==null) {
String msg="Failed to read sensor image file "+sourcePath;
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
if (debugLevel>0) System.out.println("Read "+sourcePath+" as a sensor image");
if (this.lanczos==null) generateLanczosStack();
if (this.oversampled!=sourceImageScale) generateLanczosStack(
this.lanczosA,
sourceImageScale,
this.binsPerHalfPixel);
if (impSrc.getType()!=ImagePlus.COLOR_RGB) {
String msg="Not yet implemented for this image type, only RGB24 is supported";
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
ImagePlus impOverlap=interSensor.resampleToOverlapRGB24(
true, // it is overlap of a pair of images
impSrc,
channel,
sourceImageScale,
this.lanczos,
this.binsPerHalfPixel,
maxThreads);
if (impOverlap==null){
String msg="Failed to apply overlap map to "+sourcePath;
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
if (debugLevel>1){
impOverlap.getProcessor().resetMinAndMax(); // imp_psf will be reused
impOverlap.show();
}
//resultDirectory
if (saveTiff){
String outPath=((resultDirectory.length()<1)?"":(resultDirectory+Prefs.getFileSeparator()))+impOverlap.getTitle()+".tiff";
try {
(new EyesisTiff()).saveTiff(
impOverlap,
outPath,
3, // float32
255.0,
true,
debugLevel);
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (FormatException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (ServiceException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (DependencyException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
if (convertToDouble) interSensor.argb24ToAYCbCr(
impOverlap,
iChn, // 0 - first image, 1 - second image in the pair
aYCbCr);
}
if (convertToDouble) interSensor.overlapImages=aYCbCr;
return aYCbCr;
}
public InterSensor loadPlaneMap(
String path, //
int debugLevel
){
Opener opener=new Opener();
ImagePlus impMap=opener.openImage("", path);
if (impMap==null) {
String msg="Failed to read plane-to-sensor map file "+path;
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
if (debugLevel>0) System.out.println("Read "+path+" as a plane-to-sensor map");
(new JP46_Reader_camera(false)).decodeProperiesFromInfo(impMap);
InterSensor interSensor=new InterSensor(impMap,debugLevel);
this.lastUsedInterSensor=interSensor;
for (int iChn=0;iChn1) System.out.println("loadPlaneMap():this.sensors["+channel+"]==null");
} else {
this.sensors[channel].interSensor=interSensor; // OK if some sensors are not used for the selected source files
}
}
return interSensor;
}
public void applyPlaneMap(
String path, //
String imgPathFormat,
String resultDirectory,
int sourceImageScale, // 2.0
boolean saveTiff,
int maxThreads,
int debugLevel
){
InterSensor interSensor=loadPlaneMap(
path, //
debugLevel
);
Opener opener=new Opener();
for (int iChn=0;iChn0) System.out.println("Failed to read sensor image file "+sourcePath);
continue;
}
if (debugLevel>0) System.out.println("Read "+sourcePath+" as a sensor image");
if (this.lanczos==null) generateLanczosStack();
if (this.oversampled!=sourceImageScale) generateLanczosStack(
this.lanczosA,
sourceImageScale,
this.binsPerHalfPixel);
if (impSrc.getType()!=ImagePlus.COLOR_RGB) {
String msg="Not yet implemented for this image type, only RGB24 is supported";
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
ImagePlus impOverlap=interSensor.resampleToOverlapRGB24(
false, // it is plane, not overlap of a pair of images
impSrc,
channel,
sourceImageScale,
this.lanczos,
this.binsPerHalfPixel,
maxThreads);
if (impOverlap==null){
String msg="Failed to apply plane map to "+sourcePath;
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
if (debugLevel>1){
impOverlap.getProcessor().resetMinAndMax(); // imp_psf will be reused
impOverlap.show();
}
//resultDirectory
if (saveTiff){
String outPath=((resultDirectory.length()<1)?"":(resultDirectory+Prefs.getFileSeparator()))+impOverlap.getTitle()+".tiff";
try {
(new EyesisTiff()).saveTiff(
impOverlap,
outPath,
3, // float32
255.0,
true,
debugLevel);
} catch (IOException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (FormatException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (ServiceException e) {
// TODO Auto-generated catch block
e.printStackTrace();
} catch (DependencyException e) {
// TODO Auto-generated catch block
e.printStackTrace();
}
}
}
}
public ImagePlus applyPlaneMap(
ImagePlus impSrc,
int channel,
int sourceImageScale, // 2.0
int maxThreads,
int debugLevel
){
InterSensor interSensor=this.sensors[channel].interSensor;
if (this.lanczos==null) generateLanczosStack();
if (this.oversampled!=sourceImageScale) generateLanczosStack(
this.lanczosA,
sourceImageScale,
this.binsPerHalfPixel);
if (impSrc.getType()!=ImagePlus.COLOR_RGB) {
String msg="Not yet implemented for this image type, only RGB24 is supported";
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
ImagePlus impPlane=interSensor.resampleToOverlapRGB24(
false, // it is plane, not overlap of a pair of images
impSrc,
channel,
sourceImageScale,
this.lanczos,
this.binsPerHalfPixel,
maxThreads);
if (impPlane==null){
String msg="Failed to apply plane map to channel "+channel;
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
if (debugLevel>2){
impPlane.getProcessor().resetMinAndMax(); // imp_psf will be reused
impPlane.show();
}
return impPlane;
}
public void generateAndSaveInterSensorMaps(
String directory,
String fileNameFormat, //prefix%02d-%02dsuffix
double alphaThreshold,
double overlapThreshold,
double angleSizeX,
double angleSizeY,
int overlapExtraPixels,
boolean updateStatus,
int debugLevel
){
int [][] pairs =findSensorPairs(
alphaThreshold,
overlapThreshold);
if (debugLevel>0) {
System.out.println("Detected "+pairs.length+" overlapping sensor pairs");
if (debugLevel>1){
for (int numPair=0;numPair1) System.out.println("Processing sensor pair "+(numPair+1)+ " ( of "+pairs.length+") - sensor "+
pairs[numPair][0]+" and sensor "+pairs[numPair][1]);
String title=String.format(fileNameFormat,pairs[numPair][0],pairs[numPair][1]);
ImagePlus imp=interSensorOverlapMap(
pairs[numPair][0],
pairs[numPair][1],
angleSizeX,
angleSizeY,
overlapExtraPixels,
title,
debugLevel);
if (imp!=null) {
if (debugLevel>2) {
imp.getProcessor().resetMinAndMax(); // imp_psf will be reused
imp.show();
}
String path=((directory!=null) && (directory.length()>0)?(directory+Prefs.getFileSeparator()):"")+title;
FileSaver fs=new FileSaver(imp);
String msg="Saving ovelap map for sensors "+pairs[numPair][0]+" and "+pairs[numPair][1]+": "+path;
if (updateStatus) IJ.showStatus(msg);
if (debugLevel>0) System.out.println(msg);
fs.saveAsTiffStack(path);
} else { // imp==null
if (debugLevel>0){
System.out.println("No overlap detected for sensors "+pairs[numPair][0]+" and "+pairs[numPair][1]);
}
}
}
}
/**
* Create 3xNumber of sensors (channels) map from plane pixels to individual sensor images pixels (pixel-X, pixel-Y, alpha(mask) )
* Uses intermediate equirectangular map that should be defined
* @param channels - array of sensor numbers
* @param projectionElevation - Latitude (in degrees) of the normal to the projection plane
* @param projectionYaw - Longitude (in degrees) of the normal to the projection plane
* @param projectionRoll - Rotation of the projection plane around the perpendicular from the lens centers
* @param projectionPixelSize ratio of the plane pixel size to the distance from the lens center to the projection plane
* @param projectionWidth - width of the projection rectangle
* @param projectionHeight - height of the projection rectangle
* @param projectionCenterX - X-coordinate (along the projection plane X - right) of the intersection of the projection plane with the perpendicular from the lens center
* @param projectionCenterY - Y-coordinate (along the projection plane Y - down) of the intersection of the projection plane with the perpendicular from the lens center
* @param title - Image title
* @param debugLevel - debug level
* @return image with 3*N slices - {pixelX1, pixelY1, alpha1, pixelX2, pixelY2, alpha2, ...} and metadata needed to map images
*/
// TODO: Make NaN for projectionPixelSize,projectionCenterX,projectionCenterY and 0 for projectionWidth - then find automatically
// also save as parameters
public ImagePlus getPlaneToSensorsMap(
int [] channels,
double projectionElevation,
double projectionYaw,
double projectionRoll,
double projectionPixelSize,
int projectionWidth,
int projectionHeight,
double projectionCenterX,
double projectionCenterY,
double nominalHorizontalDisparity, // nominal distance between horizontal cameras
String title,
int debugLevel
){
generateOverlapBits(0.1, channels ); // sets this.panoDegreesPerPixel among other things
// determine projectionPixelSize (should be the same for all sensors
if (Double.isNaN(projectionPixelSize) || (projectionPixelSize<=0.0)) {
projectionPixelSize=Math.PI/180.0*this.panoDegreesPerPixel;
if (debugLevel>0) System.out.println("getPlaneToSensorsMap(): Using projectionPixelSize="+projectionPixelSize);
// SensorData.EquirectangularMap erm=this.sensors[channels[0]].equirectangularMap;
// projectionPixelSize=erm.degreesPerPixel*Math.PI/180;
}
// create lat/long pair for each pixel of the projection plane
double projPsi=Math.PI/180*projectionRoll;
double projTheta=Math.PI/180*projectionElevation;
double projPhi=Math.PI/180*projectionYaw;
// rotate by - projPsi around projection plane z (normal)
double [][] aR10={
{ Math.cos(projPsi), Math.sin(projPsi), 0.0},
{-Math.sin(projPsi), Math.cos(projPsi), 0.0},
{ 0.0, 0.0, 1.0}};
Matrix R10=new Matrix(aR10,3,3);
// rotate by - projTheta around projection plane X
double [][] aR21={
{1.0, 0.0, 0.0 },
{0.0, Math.cos(projTheta), Math.sin(projTheta)},
{0.0,-Math.sin(projTheta), Math.cos(projTheta)}};
Matrix R21=new Matrix(aR21,3,3);
// rotate by -projPhi around projection plane Y
double [][] aR32={
{ Math.cos(projPhi), 0.0, Math.sin(projPhi)},
{ 0.0, 1.0, 0.0 },
{-Math.sin(projPhi), 0.0, Math.cos(projPhi)}};
Matrix MA=((new Matrix(aR32,3,3)).times(R21)).times(R10);
float [][][] projLatLong = new float [projectionHeight][projectionWidth][2];
double [][] V={{0.0},{0.0},{1.0}};
Matrix MV=new Matrix(V,3,1);
for (int py=0;py2) {
int i0=0;
double [][] dbgVect=new double [3][projectionHeight*projectionWidth];
for (int py=0;py channel2 vector
* @param channels array of channels (sensors) to use
* @param channel1 right sub-camera
* @param channel2 left sub-camera
* @param debugLevel debug level
* @param strictHorizontalDisparity - set roll to make two cameras have pure horizontal disparity (false - minimize overall roll)
* @return (Yaw, elevation, roll) of the projection plane in degrees
*/
public double [] suggestProjectionPlaneYawElevationRoll(
int [] channels,
int channel1,
int channel2,
boolean strictHorizontalDisparity,
int debugLevel
){
int debugThreshold=1;
double [][] VZ={{0.0},{0.0},{1.0}};
double [][] VX={{1.0},{0.0},{0.0}};
double [][] VY={{0.0},{1.0},{0.0}}; // up?
double [] averageViewVector={0.0,0.0,0.0};
double [] averageViewXVector={0.0,0.0,0.0};
for (int iSens=0;iSens=debugThreshold) {
System.out.println("iSens="+iSens+", chn="+chn);
System.out.println("this.sensors["+chn+"].theta"+this.sensors[chn].theta+", theta="+theta);
System.out.println("this.sensors["+chn+"].azimuth"+this.sensors[chn].azimuth+
" this.sensors["+chn+"].phi"+this.sensors[chn].phi+" (sum="+
(this.sensors[chn].azimuth+this.sensors[chn].phi)+
", phi="+phi);
}
// rotate by - projPsi around projection plane z (normal)
double [][] aR10={
{ Math.cos(psi), Math.sin(psi), 0.0},
{-Math.sin(psi), Math.cos(psi), 0.0},
{ 0.0, 0.0, 1.0}};
Matrix R10=new Matrix(aR10,3,3);
// rotate by - projTheta around projection plane X
double [][] aR21={
{1.0, 0.0, 0.0 },
{0.0, Math.cos(theta), Math.sin(theta)},
{0.0,-Math.sin(theta), Math.cos(theta)}};
Matrix R21=new Matrix(aR21,3,3);
// rotate by -projPhi around projection plane Y
double [][] aR32={
{ Math.cos(phi), 0.0, Math.sin(phi)},
{ 0.0, 1.0, 0.0 },
{-Math.sin(phi), 0.0, Math.cos(phi)}};
Matrix MA=(new Matrix(aR32,3,3)).times(R21.times(R10));
if (debugLevel>=debugThreshold) {
System.out.println("R21=");
R21.print(10,5);
System.out.println("R32=");
(new Matrix(aR32,3,3)).print(10,5);
System.out.println("MA=");
MA.print(10,5);
}
double [] sensorViwewVector= unityVector((MA.times(new Matrix(VZ))).getRowPackedCopy());
double [] sensorViwewXVector=unityVector((MA.times(new Matrix(VX))).getRowPackedCopy());
for (int i=0;i=debugThreshold) {
System.out.println("Sensor view vector=("+sensorViwewVector[0]+","+sensorViwewVector[1]+","+sensorViwewVector[2]+")");
System.out.println("Sensor view X vector=("+sensorViwewXVector[0]+","+sensorViwewXVector[1]+","+sensorViwewXVector[2]+")");
}
averageViewVector[i]+= sensorViwewVector[i];
averageViewXVector[i]+=sensorViwewXVector[i];
}
}
double [] planeZ=unityVector(averageViewVector);
double [] planeX=unityVector(crossProduct3(planeZ,crossProduct3(averageViewXVector,planeZ)));
double [] planeXM={-planeX[0],-planeX[1],-planeX[2]};
if (debugLevel>=debugThreshold) {
System.out.println("suggestProjectionPlaneYawElevationRoll(channels[],"+channel1+","+channel2+")");
System.out.println("planeZ=("+planeZ[0]+","+planeZ[1]+","+planeZ[2]+")");
System.out.println("planeX=("+planeX[0]+","+planeX[1]+","+planeX[2]+")");
System.out.println("planeXM=("+planeXM[0]+","+planeXM[1]+","+planeXM[2]+")");
}
double [] v1=lensCenterVector(channel1);
double [] v2=lensCenterVector(channel2);
double [] interVector={0.0,0.0,0.0};
for (int i=0;i=debugThreshold) {
System.out.println("interProjectionUnity=("+interProjectionUnity[0]+","+interProjectionUnity[1]+","+interProjectionUnity[2]+")");
System.out.println("horProjectionUnity=("+horProjectionUnity[0]+","+horProjectionUnity[1]+","+horProjectionUnity[2]+")");
System.out.println("vertProjectionUnity=("+vertProjectionUnity[0]+","+vertProjectionUnity[1]+","+vertProjectionUnity[2]+")");
System.out.println("cosPPRoll="+cosPPRoll+" sinPPRoll="+sinPPRoll);
System.out.println("ppYaw="+ppYaw+" ppElevation="+ppElevation+" ppRoll="+ppRoll);
}
return result;
}
/**
* Project lens center (entrance pupil) on the projection plane (going through the camera center)
* @param channel subcamera number
* @param ppYaw projection plane normal yaw in degrees, positive - CW looking from top
* @param ppElevation projection plane normal elevation in degrees (up - positive)
* @param ppRoll projection plane rotation around the normal, positive CW looking from the camera center
* @param debugLevel debug level
* @return triplet of coordinates aligned to the projection plane (z - away from the camera center)
*/
double [] projectPupil(
int channel,
double ppYaw,
double ppElevation,
double ppRoll){
double psi=Math.PI/180*ppRoll;
double theta=Math.PI/180*ppElevation; // radians
double phi=Math.PI/180*ppYaw; // radians
double [] v=lensCenterVector(channel);
double [][] V={{v[0]},{v[1]},{v[2]}};
/*1) rotate by -phi) around
| cos(phi) 0 -sin(phi) |
| 0 1 0 |
| sin(phi) 0 cos(phi) |
*/
double [][] aR0={
{Math.cos(phi),0.0,-Math.sin(phi)},
{0.0, 1.0, 0.0},
{Math.sin(phi),0.0, Math.cos(phi)}};
Matrix R0=new Matrix(aR0);
/*
2) rotate by theta around
| 1 0 0 |
| 0 cos(theta) -sin(theta) |
| 0 sin(theta) cos(theta) |
*/
double [][] aR1={
{1.0, 0.0, 0.0},
{0.0, Math.cos(theta),-Math.sin(theta)},
{0.0, Math.sin(theta), Math.cos(theta)}};
Matrix R1=new Matrix(aR1);
// 3) rotate by psi around projection plane z (normal)
double [][] aR2={
{ Math.cos(psi), -Math.sin(psi), 0.0},
{ Math.sin(psi), Math.cos(psi), 0.0},
{ 0.0, 0.0, 1.0}};
Matrix R2=new Matrix(aR2);
double [] xyz= R2.times(R1.times(R0.times(new Matrix(V)))).getRowPackedCopy();
return xyz;
}
// double [] crossProduct3(double [] a,double [] b){
/**
* Calculate sensor entrance pupil coordinates in camera coordinate system
* @param chn sensor number
* @return (x,y,z)
*/
public double [] lensCenterVector(int chn){
double theta=Math.PI/180*this.sensors[chn].theta; // radians
double phi=Math.PI/180*this.sensors[chn].phi; // radians
double azimuth=Math.PI/180*this.sensors[chn].azimuth; // radians
/* 0) Translate by distance to entrance pupil (lens center)
| Xc0 | | 0 | |Xc|
| Yc0 | = | 0 | + |Yc|
| Zc0 | | entrancePupilForward | |Zc|
*/
double [][] aT0={{0.0},{0.0},{this.sensors[chn].entrancePupilForward}};
Matrix T0=new Matrix(aT0);
/*
2) rotate by - theta around C1X:Vc2= R2*Vc1
| Xc2 | | 1 0 0 | |Xc1|
| Yc2 | = | 0 cos(theta) sin(theta) | * |Yc1|
| Zc2 | | 0 -sin(theta) cos(theta) | |Zc1|
*/
double [][] aR2={
{1.0,0.0,0.0},
{0.0,Math.cos(theta),Math.sin(theta)},
{0.0,-Math.sin(theta),Math.cos(theta)}};
Matrix R2=new Matrix(aR2);
/*
3) rotate by -(azimuth+phi) around C2Y:Vc3= R3*Vc2
| Xc3 | | cos(azimuth+phi) 0 sin(azimuth+phi) | |Xc2|
| Yc3 | = | 0 1 0 | * |Yc2|
| Zc3 | | -sin(azimuth+phi) 0 cos(azimuth+phi) | |Zc2|
*/
double [][] aR3={
{Math.cos(phi+azimuth),0.0,Math.sin(phi+azimuth)},
{0.0,1.0,0.0},
{-Math.sin(phi+azimuth),0.0,Math.cos(phi+azimuth)}};
Matrix R3=new Matrix(aR3);
/*
4) Now axes are aligned, just translate to get to eyesis coordinates: Vey= T1+Vc3
| Xey | | r * sin (azimuth) | |Xc3|
| Yey | = | height+centerAboveHorizontal | + |Yc3|
| Zey | | r * cos (azimuth) | |Zc3|
*/
double [][] aT1={
{this.sensors[chn].radius*Math.sin(azimuth)},
{this.sensors[chn].height},
{this.sensors[chn].radius*Math.cos(azimuth)}};
Matrix T1=new Matrix(aT1);
// MA=R3*R2;
// MB=T1+R3*R2*T0;
// Matrix MA=R3.times(R2);
// Matrix MB=T1.plus(R3.times(R2.times(T0)));
double [] aB=T1.plus(R3.times(R2.times(T0))).getRowPackedCopy();
return aB;
}
/**
* Create a 6-slice map of the flat rectangular intersection area of two sensors
* to the sensor pixel coordinates. The map is centered around the mid-point of the two
* sensors on the sphere, horizontal dimension (x) from sensor1 to sensor 2.
* this.panoDegreesPerPixel, this.panoLongitudeLeft, this.panoLongitudeRight,
* this.panoLatitudeTop and this.panoLatitudeBottom should be already set (i.e. by
* overlapPairsAreas() ).
* @param channel1 first sensor channel number
* @param channel2 second sensor channel number
* @param angleSizeX overlap dimension in degrees (usually narrow)
* @param angleSizeY overlap dimension in degrees (usually high)
* @param extraPixels trim result to have only this number of non-overlapping pixels along X-axis
* @param debugLevel debug level
* @return image with 6 slices - {pixelX1, pixelY1, alpha1, pixelX2, pixelY2, alpha2}
*/
public ImagePlus interSensorOverlapMap(
int channel1,
int channel2,
double angleSizeX,
double angleSizeY,
int extraPixels,
String title,
int debugLevel
){
// double [] interDistance={0.0};
double [][] vectors = new double [2][];
float [][][] latLong= mapInterSensorToLatLong(
channel1,
channel2,
angleSizeX,
angleSizeY,
vectors,
debugLevel);
int height=latLong.length;
int width=latLong[0].length;
float [][] pixels1=mapRectangleToSensor(
channel1,
latLong, //(=/-90, +/- 180)
debugLevel);
float [][] pixels2=mapRectangleToSensor(
channel2,
latLong, //(=/-90, +/- 180)
debugLevel);
// find first and last non-empty lines
// Trim result arrays
int minX=width,maxX=-1,minY=height,maxY=-1;
for (int iy=0;iy0.0) && (pixels2[2][index]>0.0)){ // overlap
if (ix>maxX) maxX=ix;
if (ixmaxY) maxY=iy;
if (iymaxX){
if (debugLevel>0){
System.out.println ("interSensorOverlapMap(): No overlap width="+width+" height="+height+
" minX="+minX+" maxX="+maxX+" minY="+minY+" maxY="+maxY);
}
return null;
}
if (debugLevel>1){
System.out.println ("interSensorOverlapMap(): width="+width+" height="+height+
" minX="+minX+" maxX="+maxX+" minY="+minY+" maxY="+maxY);
}
minX-=extraPixels;
if (minX<0) minX=0;
maxX+=extraPixels;
if (maxX>=width) maxX=width-1;
int projectionCenterX=(width-1)/2-minX; // width is odd, symmetrical around the center
int projectionCenterY=(height-1)/2-minY; // height is odd, symmetrical around the center
int oWidth=maxX-minX+1;
int oHeight=maxY-minY+1;
float [][] pixels=new float [6][oWidth*oHeight];
for (int iy=0;iy=..>+180 for longitude
* @param channel1 - number of the first camera
* @param channel1 - number of the second camera
* @param angleSizeX angular overlap size along the line connecting two cameras, in degrees
* @param angleSizeY angular overlap size perpendicular to the line connecting two cameras, in degrees
* @param vectors should be initialized as double [2][]. Returns 2 vectors:
* first - unity vector, normal to the image plane, second - vector from camera1 to camera 2, mm.
* Image plane axis X is a projection of this vector on the image plane
* @param debugLevel debug level
* @return [pixels perpendicular to line between cameras][pixels in the direction from camera 1 to camera 2]{lat, long}
*/
public float [][][] mapInterSensorToLatLong(
int channel1,
int channel2,
double angleSizeX,
double angleSizeY,
double [][] vectors, // should be initialized as double[2][]
int debugLevel
){
double lat1= this.sensors[channel1].theta; //latitude (elevation) of the first camera optical axis, degrees
double long1 = this.sensors[channel1].azimuth+this.sensors[channel1].phi; //longitude (heading+azimuth) of the first camera (CW from top) adding multiple 360 is OK here
double lat2 = this.sensors[channel2].theta; //latitude of the second camera
double long2= this.sensors[channel2].azimuth+this.sensors[channel2].phi; // longitude of the second camera, dding multiple 360 is OK here
double degreesPerPixel= this.panoDegreesPerPixel; //angular size of a pixel on equator, in degrees
// camera coordinates y - up, to the target (from the camera at zero angle/zero heading), x - right (looking to the target)
// converting sensor vector {0,0,1} to the camera coordinates. skipping roll
double theta1= lat1* Math.PI/180;
double phi1= long1*Math.PI/180;
double theta2= lat2* Math.PI/180;
double phi2= long2*Math.PI/180;
double [][] aSensorAxis={{0.0},{0.0},{1.0}};
Matrix mSensorAxis=new Matrix(aSensorAxis);
// rotate by - theta1 around sensor1 x:Vc2= R2*Vc1
double [][] aR21={
{1.0, 0.0, 0.0 },
{0.0, Math.cos(theta1), Math.sin(theta1)},
{0.0,-Math.sin(theta1), Math.cos(theta1)}};
Matrix R21=new Matrix(aR21);
// rotate by -(phi1) around C2Y:Vc3= R3*Vc2
double [][] aR31={
{ Math.cos(phi1), 0.0, Math.sin(phi1)},
{ 0.0, 1.0, 0.0 },
{-Math.sin(phi1), 0.0, Math.cos(phi1)}};
Matrix MA1=(new Matrix(aR31)).times(R21);
double [] u1=(MA1.times(mSensorAxis)).getColumnPackedCopy(); // unity vector in the direction of the first camera axis
double [][] aR22={
{1.0, 0.0, 0.0 },
{0.0, Math.cos(theta2), Math.sin(theta2)},
{0.0,-Math.sin(theta2), Math.cos(theta2)}};
Matrix R22=new Matrix(aR22);
// rotate by -(phi1) around C2Y:Vc3= R3*Vc2
double [][] aR32={
{ Math.cos(phi2), 0.0, Math.sin(phi2)},
{ 0.0, 1.0, 0.0 },
{-Math.sin(phi2), 0.0, Math.cos(phi2)}};
Matrix MA2=(new Matrix(aR32) ).times(R22);
double [] u2=(MA2.times(mSensorAxis)).getColumnPackedCopy(); // unity vector in the direction of the first camera axis
double [] rsltUZ={
u2[0]+u1[0],
u2[1]+u1[1],
u2[2]+u1[2]};
rsltUZ=unityVector(rsltUZ); // unity vector between the views of the two cameras (maybe - find center of the overlap?
if (debugLevel>1) System.out.println("mapInterSensorToLatLong(): rsltUZ={"+rsltUZ[0]+","+rsltUZ[1]+","+rsltUZ[2]+"}");
if (debugLevel>1) System.out.println("mapInterSensorToLatLong(): radius1="+this.sensors[channel1].radius+
", azimuth1="+this.sensors[channel1].azimuth+
", height1="+this.sensors[channel1].height);
if (debugLevel>1) System.out.println("mapInterSensorToLatLong(): radius2="+this.sensors[channel2].radius+
", azimuth2="+this.sensors[channel2].azimuth+
", height2="+this.sensors[channel2].height);
// lens centers locations:
double [] camera1XYZ={
this.sensors[channel1].radius*Math.sin(Math.PI/180.0* this.sensors[channel1].azimuth),
this.sensors[channel1].height,
this.sensors[channel1].radius*Math.cos(Math.PI/180.0* this.sensors[channel1].azimuth)};
double [] camera2XYZ={
this.sensors[channel2].radius*Math.sin(Math.PI/180.0* this.sensors[channel2].azimuth),
this.sensors[channel2].height,
this.sensors[channel2].radius*Math.cos(Math.PI/180.0* this.sensors[channel2].azimuth)};
if (debugLevel>1) System.out.println("mapInterSensorToLatLong(): camera1XYZ={"+camera1XYZ[0]+","+camera1XYZ[1]+","+camera1XYZ[2]+"}");
if (debugLevel>1) System.out.println("mapInterSensorToLatLong(): camera2XYZ={"+camera2XYZ[0]+","+camera2XYZ[1]+","+camera2XYZ[2]+"}");
double [] interCamerVector={
camera2XYZ[0]-camera1XYZ[0],
camera2XYZ[1]-camera1XYZ[1],
camera2XYZ[2]-camera1XYZ[2]};
if (debugLevel>1) System.out.println("mapInterSensorToLatLong(): interCamerVector={"+interCamerVector[0]+","+interCamerVector[1]+","+interCamerVector[2]+"}");
// here x,y,z - left!
double [] rsltUY=unityVector(crossProduct3(rsltUZ,interCamerVector));
// double [] rsltUX= crossProduct3(interCamerVector,rsltUY);
double [] rsltUX= crossProduct3(rsltUY,rsltUZ);
if (debugLevel>1) System.out.println("mapInterSensorToLatLong(): rsltUY={"+rsltUY[0]+","+rsltUY[1]+","+rsltUY[2]+"}");
if (debugLevel>1) System.out.println("mapInterSensorToLatLong(): rsltUX={"+rsltUX[0]+","+rsltUX[1]+","+rsltUX[2]+"}");
// calculate output dimensions
double outPixelSize=Math.PI/180.0*degreesPerPixel; // radius =1.0
double linearXHalfRange=Math.tan(angleSizeX*Math.PI/360.0)/outPixelSize; // using half-angle
double linearYHalfRange=Math.tan(angleSizeY*Math.PI/360.0)/outPixelSize; // using half-angle
int xHalfSize=(int) Math.round(linearXHalfRange);
int yHalfSize=(int) Math.round(linearYHalfRange);
float [][][] result = new float [2*yHalfSize+1][2*xHalfSize+1][2];
for (int iy=0;iy<=2*yHalfSize;iy++) for (int ix=0;ix<=2*xHalfSize;ix++){
// double y=(iy-yHalfSize)*outPixelSize;
double y=(yHalfSize-iy)*outPixelSize; // in pixel order - y-down
double x=(ix-xHalfSize)*outPixelSize;
double [] v={ // point on the result plane in camera coordinates
x*rsltUX[0]+y*rsltUY[0]+rsltUZ[0],
x*rsltUX[1]+y*rsltUY[1]+rsltUZ[1],
x*rsltUX[2]+y*rsltUY[2]+rsltUZ[2],
};
// convert to latitude/longitude
v=unityVector(v); // re-normalize to increase precision?
result[iy][ix][0]=(float) ((180.0/Math.PI)*Math.asin(v[1])); // latitude v[1] (y) - up -90..+90
result[iy][ix][1]=(float) ((180.0/Math.PI)*Math.atan2(v[0],v[2])); // longitude v[0](x) right, V[2](z) to the target +/-180
}
vectors[0]=rsltUZ;
vectors[1]=interCamerVector;
return result;
}
/**
* Calculate cross-product of 2 vectors in 3d
* @param a 3-element vector a
* @param b 3-element vector b
* @return 3-element cross product of vectors a and b
*/
double [] crossProduct3(double [] a, double [] b){
double []cp3={
a[1]*b[2]-a[2]*b[1],
a[2]*b[0]-a[0]*b[2],
a[0]*b[1]-a[1]*b[0]};
return cp3;
}
/**
* Calculate dot-product of 2 vectors
* @param a vector a
* @param b vector b
* @return dot product of two vectors
*/
double dotProduct(double [] a, double [] b){
double dp=0;
for (int i=0;i1){
System.out.println("this.sensors["+channel+"].channel].equirectangularMap.maoWOI= {x="+
woi.x+", y="+woi.y+", width="+woi.width+", height="+woi.height+"}");
System.out.println("this.panoDegreesPerPixel="+this.panoDegreesPerPixel);
System.out.println("this.panoLatitudeTop="+this.panoLatitudeTop);
System.out.println("this.panoLongitudeLeft="+this.panoLongitudeLeft);
}
double minLat=this.panoLatitudeTop-(woi.y+woi.height)*this.panoDegreesPerPixel;
double maxLat=this.panoLatitudeTop-woi.y*this.panoDegreesPerPixel;
double minLong=this.panoLongitudeLeft+woi.x*this.panoDegreesPerPixel;
double maxLong=this.panoLongitudeLeft+(woi.x+woi.width)*this.panoDegreesPerPixel;
while (minLong<-180.0) minLong+=360.0;
while (minLong>=180.0) minLong-=360.0; // -180<=long<180
while (maxLong<-180.0) maxLong+=360.0;
while (maxLong>=180.0) maxLong-=360.0; // -180<=long<180
boolean rollover= (maxLong<=minLong);
if (debugLevel>1){
System.out.println("mapRectangleToSensor("+channel+",...): minLong="+minLong+", maxLong="+maxLong+
", minLat="+minLat+", maxLat="+maxLat+", rollover="+rollover);
}
float [][] pixelCoords =new float[3][latLong.length*latLong[0].length];
for (int n=0;n1)){
System.out.println("mapRectangleToSensor("+channel+",...): lat["+iy+"]["+ix+"]="+lat+", long["+iy+"]["+ix+"]="+lng);
}
if ((latmaxLat)) continue;
if (!rollover &&((lngmaxLong))) continue;
if ( rollover &&((lngmaxLong))) continue;
if (lng1)){
System.out.println("mapRectangleToSensor("+channel+",...): apy["+iy+"]["+ix+"]="+apy+", apx["+iy+"]["+ix+"]="+apx);
}
int iapy= (int) Math.floor(apy);
int iapx= (int) Math.floor(apx);
if ((iapy>= (woi.height-1)) || (iapy<0)) continue;
if ((iapx>= (woi.width -1)) || (iapx<0)) continue;
int index00=iapx+iapy*woi.width;
if (erm.partialMap[2][index00]<=0.0) continue; // alpha is 0 - skip
int indexX0=index00+1;
int index0Y=index00+woi.width;
int indexXY=index0Y+1;
double dapy=apy-iapy;
double dapx=apx-iapx;
if ((iy==(latLong.length/2)) && (ix==(latLong[0].length/2)) && (debugLevel>1)){
System.out.println("mapRectangleToSensor("+channel+",...): dapy["+iy+"]["+ix+"]="+dapy+", dapx["+iy+"]["+ix+"]="+dapx);
}
if ((dapx>0.0) && (erm.partialMap[2][indexX0]<=0.0)) continue;
if ((dapy>0.0) && (erm.partialMap[2][index0Y]<=0.0)) continue;
if (((dapy>0.0) || (dapx>0.0)) && (erm.partialMap[2][indexXY]<=0.0)) continue;
// bi-linear interpolate
int index=iy*latLong[0].length+ix;
for (int n=0;n>24) & 0xff); // alpha channel
double r=K255*((iPixels[i]>>16) & 0xff);
double g=K255*((iPixels[i]>> 8) & 0xff);
double b=K255*( iPixels[i] & 0xff);
double y=KR*r+KG*g+KB*b;
aYCbCr[1][i]=y;
aYCbCr[2][i]=KPB*(b-y);
aYCbCr[3][i]=KPR*(r-y);
}
return aYCbCr;
}
public double [] showThresholdedDisparity(
double [][][] disparityMap,
int zeroShift, // 1/2 corrFFTSize
double threshold
){
double [] result =new double [disparityMap.length*disparityMap[0].length];
for (int index=0;indexmax){
iMax=i;
max=disparityMap[y][x][i];
}
if (max>=threshold){
if ((iMax==zeroShift) || (iMax==(disparityMap[y][x].length-1))){
result[index]=iMax-zeroShift;
} else {
// y=ax^2+bx+c, x==0 for the middle point, c=y[0] already defined
double b=0.5*(disparityMap[y][x][iMax+1]-disparityMap[y][x][iMax-1]);
double a=0.5*(disparityMap[y][x][iMax+1]+disparityMap[y][x][iMax-1])-disparityMap[y][x][iMax]; // negative
// double x=-0.5*b/a;
result[index]=iMax-zeroShift+(-0.5*b/a);
}
}
}
}
return result;
}
public double [] showThresholdedDisparityPixels(
double [][][] disparityMap,
int zeroShift, // 1/2 corrFFTSize
double threshold
){
double [] map=showThresholdedDisparity(
disparityMap,
zeroShift, // 1/2 corrFFTSize
threshold);
double [] result =new double [this.mapWidth*this.mapHeight];
for (int index=0;index0.0 when calculating areas
* @param debugLevel
* @return For each tile (with the step corrFFTSize/overlapFraction) where defined area exceeds minimal limits return horizontal profile. For undefined tiles retuern null
*/
public double [][] correlateOneRow(
DoubleFHT doubleFHT,
boolean side,
boolean autoCorrelation,
int corrFFTSize,
int overlapFraction,
int corrYC,
int maxDisparity,
double phaseCorrelationFraction,
double corrCbWeight,
double corrCrWeight,
double correlationHighPassSigma,
double correlationLowPassSigma,
double noiseNormalizationSignaY,
double noiseNormalizationSignaCbCr,
double minFracArea,
boolean useBinaryAlpha,
int debugLevel){
int numLayers=this.overlapImages.length/2;
// int indexAlpha=0;
// int indexY= 1;
int indexCb=2;
int indexCr=3;
int numTiles=this.mapWidth*overlapFraction/corrFFTSize;
if (doubleFHT==null) doubleFHT=new DoubleFHT();
double[] hamming1d=doubleFHT.getHamming1d(corrFFTSize,0.0); // pure shifted cosine, not real Hamming
double [] componentWeights={
1.0/(1.0+corrCbWeight+corrCrWeight),
corrCbWeight/(1.0+corrCbWeight+corrCrWeight),
corrCrWeight/(1.0+corrCbWeight+corrCrWeight)};
double [][][] selection = new double[numTiles] [][];
double [] corr=new double[corrFFTSize*corrFFTSize];
double [] frequencyFilter=null;
for (int nTile=0;nTile0.0) || (correlationLowPassSigma>0.0))){
frequencyFilter=doubleFHT.createFrequencyFilter(
selection[nTile][nComp+nImage*numLayers], // just for length
correlationHighPassSigma,
correlationLowPassSigma);
}
}
}
}
}
double [][] result = new double [numTiles][];
for (int nTile=0;nTile2) debugCorr=new double[span][];
for (int nTile=0;nTile2) for (int i=0;i2) for (int i=0;i=0) && (selection[numFirst][0]!=null) && (selection[numFirst][numLayers]!=null)){
for (int i=0;i0){
double rms=hFNoiseNormalize(
first, // double [] data, will be modified
sigma, // double filterSigma,
true); // boolean centerOnly);
if (debugLevel>3){
System.out.println("Correlation RMS for component "+((nComp==1)?"Y":((nComp==2)?"Cb":"Cr"))+ " was "+rms);
}
}
for (int i=0;i=0) && (dstIndex2) {
debugCorr[nDisp]=corr.clone();
System.out.println("correlateOneRow(): nTile="+nTile+" nDisp="+nDisp+" numFirst="+numFirst+" numSecond="+numSecond);
}
}
if (debugLevel>2) {
(new ShowDoubleFloatArrays()).showArrays(
debugCorr,
corrFFTSize,
corrFFTSize,
true,
"corr"+nTile+"-y"+corrYC+"-MAXDSP"+maxDisparity+"-PC"+phaseCorrelationFraction);
}
}
return result;
}
/**
* Merge several representations of the linear features (pointed from different cells)
* @param corrFFTSize size of the square tile side
* @param overlapFraction Shift by this fraction of the corrFFTSize between tiles
* @param features array of [row][column][parameters] describing the linear feature (see tileLineDetect() method)
* @param strengthMode scale normalized linear feature to (0.0): absolute strength, 1.0 - relative strength (0.01 - several merged)
* @param scaleDistances // ~1.05? compensate for decrease in measured distances to the features (partially caused by windowing)
* @param debugLevel Debug level
* @return array of the same format as input "featured" with some features merged and some - removed
*/
public double [][][] mergeLinearFeatures(
int corrFFTSize,
int overlapFraction, // default 8
double [][][] features,
double strengthMode, // 0.0 - absolute, 1.0 - relative, 0.5 - "balanced", "-1" - use phaseStrength instead
double directionTolerance,
double normalDistanceTolerance,
double tangentialDistanceTolerance,
double hostsTolerance,
double directionSigma, // make a fixed fraction of directionTolerance?
double normalDistanceSigma, // make a fixed fraction of distanceTolerance?
double tangentialDistanceSigma, // make a fixed fraction of distanceTolerance?
int minMerged,
double cellUsageShift,
double scaleDistances,
double swapTangentialTolerance,
int swapSearchRange,
int debugRow,
int debugColumn,
int debugLevel){
// int debugRow=144; //106; //153; //148;
// int debugColumn=43; //60;// 49; //53;
// int modDebugLevel=debugLevel;
int height=features.length;
int width=features[0].length;
// now include not just same cell, but all other hosts
int [][] cellUsage= getCellUsage(
corrFFTSize,
overlapFraction, // default 8
features,
directionTolerance,
normalDistanceTolerance,
tangentialDistanceTolerance,
hostsTolerance,
scaleDistances,
debugRow,
debugColumn,
debugLevel);
// find all cells where feature is in the same cell
if (debugLevel>1){
System.out.println("mergeLinearFeatures()1A:"+
" directionTolerance="+ directionTolerance+
" directionSigma="+ directionSigma+
// " kAngle="+ kAngle+
" normalDistanceTolerance="+ normalDistanceTolerance+
" tangentialDistanceTolerance="+ tangentialDistanceTolerance+
" normalDistanceSigma="+ normalDistanceSigma+
" tangentialDistanceSigma="+ tangentialDistanceSigma
// +" kDist="+ kDist
);
}
double [][][] shares=new double [height][width][];
int [][][][] hosts= new int [height][width][][];
for (int row=0;row=0) && (swapTangentialTolerance>=0)){
for (int n=1;ndebugThreshold){
System.out.println("swapHosts()1: row="+row+" col="+col+
" row0="+row0+" col0="+col0+
" dX="+dX+" dY="+dY+
" angle="+IJ.d2s(180/Math.PI*angle,2)+
" distance="+distance);
}
numOutsideCell++;
if (distance<0.0){ // maybe it is not needed here at all?
distance=-distance;
angle+=Math.PI;
angle-=2*Math.PI*Math.round(angle/(2*Math.PI));
}
if ((distance>maximalDistance) && // just for debug/statistics, remove border elements
(row> overlapFraction/2) &&
(row< (height-overlapFraction/2)) &&
(col> overlapFraction/2) &&
(col< (width-overlapFraction/2))) {
maximalDistance=distance;
if (modDebugLevel>debugThreshold){
System.out.println("swapHosts()2: row="+row+" col="+col+ " distance="+distance+" maximalDistance"+maximalDistance);
}
}
double sinA=Math.sin(angle);
double cosA=Math.cos(angle);
for (int dRow=-scanRange;dRow<=scanRange;dRow++) {
int row1=row0+dRow;
if ((row1>=0) && (row1=0) && (col1debugThreshold){
System.out.println("swapHosts()3: row="+row+" col="+col+
" row0="+row0+" col0="+col0+
" row1="+row1+" col1="+col1+
" cell is "+((features[row1][col1]== null)?"FREE":"USED")+
" xc="+xc+" yc="+yc+
" xc*cosA+yc*sinA="+(xc*cosA+yc*sinA)+
" row1="+row1+" col1="+col1+
" angle="+IJ.d2s(180/Math.PI*angle,2)+
" distance="+distance+
" newDist="+newDist);
}
}
}
}
}
}
}
if (bestNewDisthostStrengths[bestColRow[1]][bestColRow[0]]){
hostStrengths[bestColRow[1]][bestColRow[0]]=thisStrength;
hosts[bestColRow[1]][bestColRow[0]][0]=col;
hosts[bestColRow[1]][bestColRow[0]][1]=row;
}
if (modDebugLevel>debugThreshold){
System.out.print ("swapHosts()4: bestColRow[1]="+bestColRow[1]+" bestColRow[0]="+bestColRow[0]);
System.out.print (
" hosts["+bestColRow[1]+"]["+bestColRow[0]+"][0]="+hosts[bestColRow[1]][bestColRow[0]][0]+
" hosts["+bestColRow[1]+"]["+bestColRow[0]+"][1]="+hosts[bestColRow[1]][bestColRow[0]][1]);
System.out.println(" thisStrength="+thisStrength+" hostStrengths[bestColRow[1]][bestColRow[0]]="+hostStrengths[bestColRow[1]][bestColRow[0]]);
}
}
}
}
}
if (debugLevel>2) System.out.println(">>>> swapHosts(): numWantMove="+numWantMove+" numOutsideCell="+numOutsideCell+" maximaDistance="+maximalDistance);
int [] result={numWantMove,numWantMove,numOutsideCell};
if (numWantMove==0) return result;
// now move linear features to closer to them cells
int numMoved=0;
for (int row=0;rowdebugThreshold){
int corrYC= row*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC= col*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrYC1=row1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC1=col1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.print("swapHosts()5: "+" YC="+corrYC+" XC="+corrXC+
" row="+row+" col="+col+
" row1="+row1+" col1="+col1+" YC1="+corrYC1+" XC1="+corrXC1+
" angle="+IJ.d2s(180/Math.PI*angle,2)+
" distance="+distance);
}
// distance=xc*cosA+yc*sinA - distance;
distance-=xc*cosA+yc*sinA;
if (distance<0.0){ // maybe it is not needed here at all?
distance=-distance;
angle+=Math.PI;
angle-=2*Math.PI*Math.round(angle/(2*Math.PI));
feature[5]=-feature[5];
}
feature[3]=angle;
feature[4]=distance;
if (modDebugLevel>debugThreshold){
System.out.println( " new angle="+IJ.d2s(180/Math.PI*angle,2)+
" new distance="+distance);
}
features[row][col]=feature;
cellUsage[row][col]=cellUsage[row1][col1];
cellUsage[row1][col1]=-1;
numMoved++;
}
}
result[0]=numMoved;
if (debugLevel>1) System.out.println(">>>> swapHosts(): numMoved="+numMoved+" numWantMove="+numWantMove+" numOutsideCell="+numOutsideCell+" maximaDistance="+maximalDistance);
return result;
}
private double getFeatureMatch(
int row10,
int col10,
int row20,
int col20,
int gridStep,
double [][][] features,
double directionTolerance,
double normalDistanceTolerance,
double tangentialDistanceTolerance,
double directionSigma, // make a fixed fraction of directionTolerance?
double normalDistanceSigma, // make a fixed fraction of distanceTolerance?
double tangentialDistanceSigma, // make a fixed fraction of distanceTolerance?
double scaleDistances,
int debugLevel){
// swap ends to make sure the result absolutely does not depend on the order in a pair
int row1,row2,col1,col2;
if ((row20>row10) || ((row20==row10) && (col20> col10))){
row1=row10;
col1=col10;
row2=row20;
col2=col20;
} else {
row1=row20;
col1=col20;
row2=row10;
col2=col10;
}
int height=features.length;
int width= features[0].length;
if ( (row1<0) || (col1<0) || (row1>=height) || (col1>=width) ||
(row2<0) || (col2<0) || (row2>=height) || (col2>=width) ||
(features[row1][col1]==null) || (features[row2][col2]==null)) return 0;
double angle1= features[row1][col1][3];
angle1-=(2*Math.PI)*Math.round(angle1/(2*Math.PI));
double angle2= features[row2][col2][3];
angle2-=(2*Math.PI)*Math.round(angle2/(2*Math.PI));
double dAngle=angle2-angle1;
dAngle-=Math.PI*Math.round(dAngle/Math.PI);
if (Math.abs(dAngle)>directionTolerance) return 0;
if (Math.abs(angle1-angle1)>(Math.PI/2)){
angle2+=Math.PI;
angle2-=(2*Math.PI)*Math.round(angle2/(2*Math.PI));
}
double angle=(angle1+angle2)/2; // here ignoring weights
double sinAngle=Math.sin(angle);
double cosAngle=Math.cos(angle);
double distance1= scaleDistances*features[row1][col1][4];
double distance2= scaleDistances*features[row2][col2][4];
double x1=distance1*Math.cos(angle1);
double y1=distance1*Math.sin(angle1);
double x2=distance2*Math.cos(angle2)+gridStep*(col2-col1);
double y2=distance2*Math.sin(angle2)+gridStep*(row2-row1);
double dX=x2-x1;
double dY=y2-y1;
double dNormal=-sinAngle*dX+cosAngle*dY;
if (Math.abs(dNormal)>normalDistanceTolerance) return 0;
double dTangential=cosAngle*dX+sinAngle*dY;
if (Math.abs(dTangential)>tangentialDistanceTolerance) return 0;
if ((directionSigma<=0) || (normalDistanceSigma<=0) || (tangentialDistanceSigma<=0)) return 1.0;
return Math.exp(-(
dAngle*dAngle/(directionSigma*directionSigma)+
dNormal*dNormal/(normalDistanceSigma*normalDistanceSigma)+
dTangential*dTangential/(tangentialDistanceSigma*tangentialDistanceSigma)
)/2.0);
}
private double [][][] combineFeatures(
int [][][][] hosts,
double [][][] shares,
int [][] cellUsage,
int corrFFTSize,
int overlapFraction, // default 8
double [][][] features,
double strengthMode, // 0.0 - absolute, 1.0 - relative, 0.5 - "balanced", "-1" - use phaseStrength instead
double directionTolerance,
double normalDistanceTolerance,
double tangentialDistanceTolerance,
double hostsTolerance,
double directionSigma, // make a fixed fraction of directionTolerance?
double normalDistanceSigma, // make a fixed fraction of distanceTolerance?
double tangentialDistanceSigma, // make a fixed fraction of distanceTolerance?
double scaleDistances,
int debugRow,
int debugColumn,
int debugLevel){
int modDebugLevel=debugLevel;
int height=features.length;
int width=features[0].length;
int scanAroundCells=overlapFraction/2;
int gridStep=corrFFTSize/overlapFraction;
double [][][] filteredFeatures = new double [height][width][];
double interCenter=corrFFTSize/overlapFraction;
// find all cells where feature is in closer to the cell center, than matching feature - to the center of any other cell
// int [][] theseHosts= new int [(2*scanAroundCells+1)*(2*scanAroundCells+1)][2];
// double[] theseShares=new double [(2*scanAroundCells+1)*(2*scanAroundCells+1)];
for (int row=0;row0){ // "hosts" only
// modDebugLevel=debugLevel+(((row==debugRow) && (col==debugColumn))?2:0);
modDebugLevel=debugLevel+(((Math.abs(row-debugRow)<=1) && (Math.abs(col-debugColumn)<=1))?2:0);
double distance=scaleDistances*features[row][col][4];
double angle= features[row][col][3];
if (distance<0.0){ // maybe it is not needed here at all?
distance=-distance;
angle+=Math.PI;
angle-=2*Math.PI*Math.round(angle/(2*Math.PI));
// phaseAtZero=-phaseAtZero; // not used here
}
double dX=distance*Math.cos(angle);
double dY=distance*Math.sin(angle);
int row0=row+ ((int) Math.round(dY/gridStep));
int col0=col+ ((int) Math.round(dX/gridStep));
if (modDebugLevel>2){
int corrYC=row*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println("combineFeatures()0:, row="+row+" col="+col+" YC="+corrYC+" XC="+corrXC+
" row0="+row0+" col0="+col0+
" angle="+IJ.d2s(180/Math.PI*angle,2)+
" distance="+IJ.d2s(distance,2)+
" phase zero="+IJ.d2s(180/Math.PI*features[row][col][5],2)
);
}
double sumWeights=0.0, sumWDist=0.0, sumX=0.0, sumY=0.0, sumSinAngle=0.0, sumCosAngle=0.0, sumSinPhase=0.0, sumCosPhase=0.0,
sumAbsStrength=0.0, sumReferenceLevel=0.0, sumPhaseStrength=0.0;
for (int dRow=-scanAroundCells;dRow<=scanAroundCells;dRow++) {
int row1=row0+dRow;
if ((row1>=0) && (row1=0) &&(col1=0) share=shares[row1][col1][neibIndex]; // so the "energy" will be split between the hosts, preserving totals
} else if ((row==row1) && (col==col1)){
share=1.0; // self
}
if (share>0.0){
double wDistance=getFeatureMatch(
row,
col,
row1,
col1,
gridStep,
features,
directionTolerance,
normalDistanceTolerance,
tangentialDistanceTolerance,
directionSigma, // make a fixed fraction of directionTolerance?
normalDistanceSigma, // make a fixed fraction of distanceTolerance?
tangentialDistanceSigma, // make a fixed fraction of distanceTolerance?
scaleDistances,
debugLevel);
if (wDistance>0){
double [] feature=features[row1][col1].clone();
// compensate for negative distance
if (feature[4]<0.0){
feature[4]=-feature[4];
feature[3]+=Math.PI;
feature[3]-=2*Math.PI*Math.round(feature[3]/(2*Math.PI));
feature[5]=-feature[5]; // not used here
}
double dX1=scaleDistances*feature[4]*Math.cos(feature[3])+(interCenter* (col1-col));// dCol);
double dY1=scaleDistances*feature[4]*Math.sin(feature[3])+(interCenter* (row1-row));// dRow);
double dAngle=feature[3]-angle;
dAngle-=2*Math.PI*Math.round(dAngle/(2*Math.PI));
if (Math.abs(dAngle)>(Math.PI/2)){
feature[3]+=Math.PI;
feature[3]-=2*Math.PI*Math.round(feature[3]/(2*Math.PI));
feature[5]=-feature[5]; //phaseAtZero
}
double wSharedDistance=wDistance * share;
double w=combinedStrength(strengthMode,features[row][col])*wSharedDistance;
if (modDebugLevel>2){
int corrYC=row1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println("combineFeatures()1:, row1="+row1+" col1="+col1+" YC="+corrYC+" XC="+corrXC+
" strength="+features[row1][col1][0]+
" refLevel="+features[row1][col1][1]+
" angle="+IJ.d2s(180/Math.PI*features[row1][col1][3],2)+
" distance="+IJ.d2s(scaleDistances*features[row1][col1][4],2)+
" phase zero="+IJ.d2s(180/Math.PI*features[row1][col1][5],2)
);
System.out.println("combineFeatures()2:, row1="+row1+" col1="+col1+" YC="+corrYC+" XC="+corrXC+
" strength="+feature[0]+
" angle="+IJ.d2s(180/Math.PI*feature[3],2)+
" distance="+IJ.d2s(scaleDistances*feature[4],2)+
" phase zero="+IJ.d2s(180/Math.PI*feature[5],2)
);
System.out.println("combineFeatures()3:, w="+w+
" wDistance="+wDistance+
" wSharedDistance="+wSharedDistance+
" dX1="+dX1+
" dY1="+dY1+
" dX="+ dX+
" dY="+ dY+
" Math.sin(feature[3])="+Math.sin(feature[3])+
" Math.cos(feature[3])="+Math.cos(feature[3])+
" Math.sin(feature[5])="+Math.sin(feature[5])+
" Math.cos(feature[5])="+Math.cos(feature[5])
);
}
sumWeights+=w;
sumX+=w*dX1;
sumY+=w*dY1;
sumSinAngle+=w*Math.sin(feature[3]); // direction
sumCosAngle+=w*Math.cos(feature[3]); // direction
sumSinPhase+=w*Math.sin(feature[5]); // phaseAtZero
sumCosPhase+=w*Math.cos(feature[5]); // phaseAtZero
sumWDist+= wSharedDistance;
sumAbsStrength+= wSharedDistance*feature[0];
sumReferenceLevel+= wSharedDistance*feature[1];
sumPhaseStrength+= wSharedDistance*feature[2];
cellUsage[row][col]++;
if (modDebugLevel>2){
System.out.println("combineFeatures()4:"+
" sumSinAngle="+sumSinAngle+
" sumCosAngle="+sumCosAngle+
" sumSinPhase="+sumSinPhase+
" sumCosPhase="+sumCosPhase+
" cellUsage["+row+"]["+col+"]="+cellUsage[row][col]);
}
} else {
int corrYC=row1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println("*** BUG: row1="+row1+" col1="+col1+" YC="+corrYC+" XC="+corrXC);
//bug
}
}
}
}
}
}
cellUsage[row][col]--; // host was counted twice
angle= Math.atan2(sumSinAngle,sumCosAngle);
double phaseAtZero=Math.atan2(sumSinPhase,sumCosPhase);
double x=sumX/sumWeights;
double y=sumY/sumWeights;
// double distance=x*Math.sin(angle)-y*Math.cos(angle);
distance=x*Math.cos(angle)+y*Math.sin(angle);
if (modDebugLevel>2){
System.out.println("combineFeatures()5: "+
" angle="+IJ.d2s(180/Math.PI*angle,2)+
" distance="+distance+
" phaseAtZero="+phaseAtZero+
" x="+x+
" y="+y);
}
if (distance<0.0){
distance=-distance;
angle+=Math.PI;
angle-=2*Math.PI*Math.round(angle/(2*Math.PI));
phaseAtZero=-phaseAtZero;
}
double [] feature={
sumAbsStrength/sumWDist,
sumReferenceLevel/sumWDist,
sumPhaseStrength/sumWDist,
angle,
distance, // here - already scaled
phaseAtZero,
cellUsage[row][col]
};
if (modDebugLevel>2){
System.out.println("combineFeatures()6: "+
" angle="+IJ.d2s(180/Math.PI*angle,2)+
" distance="+distance+
" phaseAtZero="+phaseAtZero+
" x="+x+
" y="+y);
System.out.println("combineFeatures()7:"+
" strength="+feature[0]+
" refLevel="+feature[1]+
" angle="+IJ.d2s(180/Math.PI*feature[3],2)+
" distance="+IJ.d2s(feature[4],2)+
" phase zero="+IJ.d2s(180/Math.PI*feature[5],2)
);
}
filteredFeatures[row][col]=feature;
}
}
return filteredFeatures;
}
public double [] displayUsage(
double [][][] features,
int corrFFTSize,
int overlapFraction, // default 8
double scale
){
double [] usage=new double[this.mapWidth*this.mapHeight];
int k=corrFFTSize/overlapFraction;
for (int index=0;index6) d=features[y][x][6];
else d=1.0;
}
usage[index]=d*scale;
}
return usage;
}
private void getFeatureShares(
int [][][][] hosts,
double [][][] shares,
int [][] cellUsage,
int corrFFTSize,
int overlapFraction, // default 8
double [][][] features,
double strengthMode, // 0.0 - absolute, 1.0 - relative, 0.5 - "balanced", "-1" - use phaseStrength instead
double directionTolerance,
double normalDistanceTolerance,
double tangentialDistanceTolerance,
double hostsTolerance,
double directionSigma, // make a fixed fraction of directionTolerance?
double normalDistanceSigma, // make a fixed fraction of distanceTolerance?
double tangentialDistanceSigma, // make a fixed fraction of distanceTolerance?
double scaleDistances,
int debugRow,
int debugColumn,
int debugLevel){
int modDebugLevel=debugLevel;
int height=features.length;
int width=features[0].length;
int scanAroundCells=overlapFraction/2;
int gridStep=corrFFTSize/overlapFraction;
// find all cells where feature is in closer to the cell center, than matching feature - to the center of any other cell
int [][] theseHosts= new int [(2*scanAroundCells+1)*(2*scanAroundCells+1)][2];
double[] theseShares=new double [(2*scanAroundCells+1)*(2*scanAroundCells+1)];
for (int row=0;row2){
int corrYC=row*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println(">>>>++++ getFeatureShares()1: row="+row+" col="+col+" YC="+corrYC+" XC="+corrXC+
" row0="+row0+" col0="+col0+
" distance="+distance+
" angle="+IJ.d2s(180/Math.PI*angle,2)+
" dX="+dX+
" dY="+dY);
}
for (int dRow=-scanAroundCells;dRow<=scanAroundCells;dRow++) {
int row1=row0+dRow;
if ((row1>=0) && (row1=0) && (col13){
int corrYC=row1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println("getFeatureShares()1a: row1="+row1+" col1="+col1+" YC="+corrYC+" XC="+corrXC+
" cellUsage[row1][col1]="+cellUsage[row1][col1]);
}
if (cellUsage[row1][col1]>0){ //cellUsage[row][col] - just to be sure, not needed
double wDist=getFeatureMatch(
row,
col,
row1,
col1,
gridStep,
features,
directionTolerance,
normalDistanceTolerance,
tangentialDistanceTolerance,
directionSigma, // make a fixed fraction of directionTolerance?
normalDistanceSigma, // make a fixed fraction of distanceTolerance?
tangentialDistanceSigma, // make a fixed fraction of distanceTolerance?
scaleDistances,
debugLevel);
if (wDist>0){
theseHosts[numHost][0]=col1;
theseHosts[numHost][1]=row1;
double w=combinedStrength(strengthMode,features[row1][col1])*wDist;
theseShares[numHost]=w;
sumShares+=w;
numHost++;
if (modDebugLevel>2){
int corrYC=row1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println("getFeatureShares()2: row1="+row1+" col1="+col1+" YC="+corrYC+" XC="+corrXC+
" wDist="+wDist+
" w="+w+
" sumShares="+sumShares+
" numHost="+numHost);
}
} else {
if (modDebugLevel>2){
int corrYC=row1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println("getFeatureShares()2: row1="+row1+" col1="+col1+" YC="+corrYC+" XC="+corrXC+
" wDist="+wDist+
" sumShares="+sumShares+
" numHost="+numHost);
}
}
}
}
}
}
}
if (numHost>0) {
hosts[row][col]=new int [numHost][2];
shares[row][col]=new double [numHost];
for (int i=0;i0 - number of cells combined, <0: 1+y*width+x of the "host"
for (int row=0;row2){
int corrYC=row*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println(">>>> getCellUsage(): row="+row+" col="+col+" YC="+corrYC+" XC="+corrXC+
" row0="+row0+" col0="+col0+
" dX="+IJ.d2s(dX,3)+" dY="+IJ.d2s(dY,3)+
" range="+range+
" angle="+IJ.d2s(180/Math.PI*angle,2)+
" distance="+IJ.d2s(distance,3)+
" strength="+IJ.d2s(features[row][col][0],4)+
" refLev="+IJ.d2s(features[row][col][1],4));
}
for (int dRow=-range;(dRow<=range) && best;dRow++) {
int row1=row0+dRow;
if ((row1>=0) && (row1=0) && (col10) &&(modDebugLevel>2)){
int corrYC=row1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
int corrXC=col1*corrFFTSize/overlapFraction+corrFFTSize/overlapFraction/2;
System.out.println("getCellUsage(): row1="+row1+" col1="+col1+" YC="+corrYC+" XC="+corrXC+
" w="+w+
" best="+best+
" distance="+IJ.d2s(distance,3)+
" scaleDistances*features["+row1+"]["+col1+"]["+4+"]="+scaleDistances*features[row1][col1][4]);
}
if (!best) break;
}
}
}
}
if (best) {
cellUsage[row][col]=1;
numHostCells++;
}
}
}
if (debugLevel>1){
System.out.println("mergeLinearFeatures(): Detected "+numHostCells+" host cells (common feature is closer to its center than to any other cell.");
}
return cellUsage;
}
/**
*
* @param corrFFTSize
* @param overlapFraction
* @param features
* @param ignorePhase
* @param strengthMode
* @param phaseIntegrationWidth
* @param resultHighPass
* @param threadsMax
* @param debugLevel
* @return
*/
public double [] reconstructImageFeature(
final int corrFFTSize,
final int overlapFraction, // default 8
final double [][][] features,
final boolean ignorePhase,
final boolean preserveDC,
final double strengthMode, // 0.0 - absolute, 1.0 - relative, 0.5 - "balanced", "-1" - use phaseStrength instead
final double phaseIntegrationWidth, // use the same as during extraction?
final double resultHighPass, //cutting "ribbon" near zero frequency - influences length of the detected line
final int debugRow,
final int debugColumn,
final int threadsMax,
final int debugLevel){
// final int numTileRows= this.mapHeight*overlapFraction/corrFFTSize+(((this.mapHeight*overlapFraction/corrFFTSize)*corrFFTSize/overlapFraction1)) {
// System.out.println("reconstructImageFeature(): row="+row+" col="+col+" corrYC="+corrYC+" corrXC="+corrXC);
// }
/*
if ((Math.abs(row-debugRow)<6) && (Math.abs(col-debugColumn)<6) && (feature!=null) && (debugLevel>1)){
System.out.println ("reconstructImageFeature(): tile="+tile+" ( of "+numTiles+"), row="+row+" col="+col+" corrYC="+corrYC+" corrXC="+corrXC+
" strength="+feature[0]+
" angle="+IJ.d2s(180/Math.PI*feature[3],2)+
" distance="+IJ.d2s(feature[4],2)+
" phase zero="+IJ.d2s(180/Math.PI*feature[5],2)
);
}
*/
synchronized(this){
int indexSrc=0;
for (int y=corrYC-hSize;y<(corrYC+hSize);y++){
if ((y>=0) && (y=0) && (xfeaturesPlot.length) || (indexSrc>tilePlot.length)) {
// System.out.println("reconstructImageFeature(): x="+x+" y="+y+" row="+row+" col="+col+" corrYC="+corrYC+" corrXC="+corrXC+
// " indexDst="+indexDst+" featuresPlot.length="+featuresPlot.length+" indexSrc="+indexSrc+" tilePlot.length="+tilePlot.length);
// }
featuresPlot[indexDst]+=tilePlot[indexSrc];
}
indexSrc++;
indexDst++;
}
} else indexSrc+=corrFFTSize;
}
}
final int numFinished=tilesFinished.getAndIncrement();
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
IJ.showProgress(numFinished,numTiles);
}
});
}
}
}
};
}
startAndJoin(threads);
IJ.showProgress(1.0);
if (debugLevel>1){
System.out.println("=== reconstructImageFeature(): number of non-empty cells="+numNonEmpty.get()+
" (of "+(numTileRows*numTileColumns)+"), "+100.0*numNonEmpty.get()/((numTileRows*numTileColumns))+"%");
}
return featuresPlot;
}
/**
* Plot the extracted linear feature to a square tile
* @param doubleFHT - instance of doubleFHT. if not null, will be reused and the cached data will not be re-calculated
* @param corrFFTSize size of the square tile side
* @param feature array of parameters describing the linear feature (see tileLineDetect() method)
* @param ignorePhase always plot feature with zero phase (white on black)
* @param preserveDC do not remove DC component from the plotted feature
* @param strengthMode scale normalized linear feature to (0.0): absolute strength, 1.0 - relative strength (0.00.0)?( absoluteStrength/Math.pow(referenceLevel,strengthMode)):phaseStrength);
if (ignorePhase) phaseAtZero=0.0;
if (doubleFHT==null )doubleFHT= new DoubleFHT();
int length=corrFFTSize*corrFFTSize;
int [] bandIndices= doubleFHT.getBandIndices(
preserveDC,
phaseIntegrationWidth,
angle);
double [] bandMask= doubleFHT.getRibbonMask(
preserveDC,
bandIndices,
phaseIntegrationWidth,
resultHighPass,
angle);
double [] phaseStepCorr=doubleFHT.compensatePhaseStep(
angle,
phaseAtZero,
distance,
0, //phaseTolerance, // if >0 will zero amplitude if the phase is too off - - only used with non-zero modifiedAmpPhase
null, // modifiedAmpPhase, // modifies phase and/or amplitude (if phaseTolerance>0),
bandIndices,
0); // zeroBinHalfSize); - only used with non-zero modifiedAmpPhase
double[] hamming1d=doubleFHT.getHamming1d(corrFFTSize,0.0); // complete zero out on the border
double [] uniform=new double[length];
for (int iy=0;iy2){
double [] debugBandMask=new double [corrFFTSize*corrFFTSize];
for (int i=0;i2){
(new ShowDoubleFloatArrays()).showArrays(
reconstructedFeature,
corrFFTSize,
corrFFTSize,
"reconstructedFeature");
}
return reconstructedFeature;
}
private double combinedStrength(
double strengthMode,
double [] feature){
return (strengthMode<-1.0)?1.0:((strengthMode>0.0)?( feature[0]/Math.pow(feature[1],strengthMode)):feature[2]);
}
/**
* Calculate parameters of the line features in the current image
* @param side 0- left, 1 - right image
* @param corrFFTSize size of the square tile side
* @param overlapFraction overlap tiles by this fraction of the tile size (corrFFTSize)
* @param alphaThreshold Minimal required overlap of the tile with image alpha (0 - do not calculate)
* @param useBinaryAlpha consider all alpha>0.0 as 1.0
* @param correlationHighPassSigma High pass filter before feature extraction (measured in frequency samples)
* @param correlationLowPassSigma Low-pass filter for correlation (as fraction of the frequency range)
* @param phaseIntegrationWidth - Width of the "band" when looking for the linear phase (actual will be smaller because of Hammimng window)
* @param resultHighPass (used only for strength calculation) - cut frequencies near zero for the linear phase band (shortens the result line);
* @param dispertionCost - >0 - use phase dispersion when calculating quality of phase match, 0 - only weights.
* @param featureFilter Bitmask enabling different phase half-steps (+1 - 0, +2 - pi/2, +4 - pi, +8 - 3pi/2. Value zero allows arbitrary step
* step 0 corresponds to thin white line, pi - thin black line, +pi/2 - edge black-to-white in the direction of directionAngle,
* 3pi/2 - white-to-black in the direction of directionAngle
* @param minRMSFrac Minimal frequency sample value relative to RMS to be considered when looking for the linear phase feature (0.0 - skip test)
* @param minAbs Minimal frequency sample absolute value to be considered when looking for the linear phase feature (0.0 - skip test)
* @param maxPhaseMismatch calculate phase differences for each full (above threshold) set of 4 samples (in a 2x2 square) and measure the distance
* between two centers of diagonals, discard set if this distance is above threshold (0.0 - skip test)
* @param calculateStrength If true, calculate feature strength by correlating the source image with the simulated linear phase band, if false (faster processing)
* these two fields will be returned equal 0.0
* @param threadsMax maximal number of parallel threads to use
* @param debugLevel Debug level, turning on more text/graphic output when higher
* @return array [tileY][tileX] {feature parameters (or null if feature failed to be detected)}:
* 0 - absolute strength of the feature (proportional to the image contrast)
* 1 - Strength reference level (RMS of the frequency samples amplitudes (in the selected band in the frequency domain)
* 2 - phaseStrength - composite value calculated when selecting the feature direction - this value is calculated even if calculateStrength is false
* 3 - angle from the tile center to the detected feature (perpendicular to the feature direction itself), zero - left, pi/2 - down (pixel coordiantes)
* 4 - distance in pixels from the tile center to the feature
* 5 - phaseAtZero Detected feature type (may be constrained by featureFilter - see that parameter for meaning of the different phase values
*/
public double [][][] detectLinearFeatures(
final boolean side,
final int corrFFTSize,
final int overlapFraction, // default 8
final double alphaThreshold, // minFracArea
final boolean useBinaryAlpha,
final double correlationHighPassSigma,
final double correlationLowPassSigma,
final double phaseIntegrationWidth,
final double resultHighPass, //cutting "ribbon" near zero frequency - influences length of the detected line (used only for strength calculation)
final double dispertionCost,
final int featureFilter,
final double minRMSFrac,
final double minAbs,
final double maxPhaseMismatch,
final boolean calculateStrength,
final int debugRow,
final int debugColumn,
final int threadsMax,
final int debugLevel) {
return detectLinearFeatures(
side?1:0,
corrFFTSize,
overlapFraction, // default 8
alphaThreshold, // minFracArea
useBinaryAlpha,
correlationHighPassSigma,
correlationLowPassSigma,
phaseIntegrationWidth,
resultHighPass, //cutting "ribbon" near zero frequency - influences length of the detected line (used only for strength calculation)
dispertionCost,
featureFilter,
minRMSFrac,
minAbs,
maxPhaseMismatch,
calculateStrength,
debugRow,
debugColumn,
threadsMax,
debugLevel);
}
public double [][][] detectLinearFeatures(
final int imageNumber, //boolean side,
final int corrFFTSize,
final int overlapFraction, // default 8
final double alphaThreshold, // minFracArea
final boolean useBinaryAlpha,
final double correlationHighPassSigma,
final double correlationLowPassSigma,
final double phaseIntegrationWidth,
final double resultHighPass, //cutting "ribbon" near zero frequency - influences length of the detected line (used only for strength calculation)
final double dispertionCost,
final int featureFilter,
final double minRMSFrac,
final double minAbs,
final double maxPhaseMismatch,
final boolean calculateStrength,
final int debugRow,
final int debugColumn,
final int threadsMax,
final int debugLevel){
if (debugLevel>1) System.out.println ("detectLinearFeatures("+imageNumber+",...)");
final int numTileRows= this.mapHeight*overlapFraction/corrFFTSize+(((this.mapHeight*overlapFraction/corrFFTSize)*corrFFTSize/overlapFraction2) System.out.println ("detectLinearFeatures(): tile="+tile+" ( of "+numTiles+"), row="+row+" col="+col+" corrYC="+corrYC+" corrXC="+corrXC);
result[row][col]= tileLineDetect(
doubleFHT,
imageNumber, //side,
alphaThreshold,
useBinaryAlpha,
corrFFTSize,
corrXC,
corrYC,
correlationHighPassSigma,
correlationLowPassSigma,
phaseIntegrationWidth,
resultHighPass, //cutting "ribbon" near zero frequency - influences length of the detected line (used only for strength calculation)
dispertionCost,
featureFilter,
minRMSFrac,
minAbs,
maxPhaseMismatch,
calculateStrength,
debugLevel);
final int numFinished=tilesFinished.getAndIncrement();
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
IJ.showProgress(numFinished,numTiles);
}
});
// if ((result[row][col]!=null) && (debugLevel>1)) System.out.println ("detectLinearFeatures(): tile="+tile+" ( of "+numTiles+"), row="+row+" col="+col+" corrYC="+corrYC+" corrXC="+corrXC+
// " strength="+result[row][col][0]);
if ((row==debugRow) && ((col==debugColumn)) && (result[row][col]!=null) && (debugLevel>1)){
System.out.println ("detectLinearFeatures(): tile="+tile+" ( of "+numTiles+"), row="+row+" col="+col+" corrYC="+corrYC+" corrXC="+corrXC+
" strength="+result[row][col][0]+
" angle="+IJ.d2s(180/Math.PI*result[row][col][3],2)+
" distance="+IJ.d2s(result[row][col][4],2)+
" phase zero="+IJ.d2s(180/Math.PI*result[row][col][5],2)+
" featureFilter="+featureFilter
);
}
}
}
};
}
startAndJoin(threads);
IJ.showProgress(1.0);
return result;
}
/**
* Calculate parameters of the line feature in the current image tile
* @param doubleFHT - instance of doubleFHT. if not null, will be reused and the cached data will not be re-calculated
* @param side 0- left, 1 - right image
* @param alphaThreshold Minimal required overlap of the tile with image alpha (0 - do not calculate)
* @param useBinaryAlpha consider all alpha>0.0 as 1.0
* @param corrFFTSize size of the square tile side
* @param corrXC Tile center X
* @param corrYC Tile center Y
* @param correlationHighPassSigma High pass filter before feature extraction (measured in frequency samples)
* @param correlationLowPassSigma Low-pass filter for correlation (as fraction of the frequency range)
* @param phaseIntegrationWidth - Width of the "band" when looking for the linear phase (actual will be smaller because of Hammimng window)
* @param resultHighPass (used only for strength calculation) - cut frequencies near zero for the linear phase band (shortens the result line);
* @param dispertionCost - >0 - use phase dispersion when calculating quality of phase match, 0 - only weights.
* @param featureFilter Bitmask enabling different phase half-steps (+1 - 0, +2 - pi/2, +4 - pi, +8 - 3pi/2. Value zero allows arbitrary step
* step 0 corresponds to thin white line, pi - thin black line, +pi/2 - edge black-to-white in the direction of directionAngle,
* 3pi/2 - white-to-black in the direction of directionAngle
* @param minRMSFrac Minimal frequency sample value relative to RMS to be considered when looking for the linear phase feature (0.0 - skip test)
* @param minAbs Minimal frequency sample absolute value to be considered when looking for the linear phase feature (0.0 - skip test)
* @param maxPhaseMismatch calculate phase differences for each full (above threshold) set of 4 samples (in a 2x2 square) and measure the distance
* between two centers of diagonals, discard set if this distance is above threshold (0.0 - skip test)
* @param calculateStrength If true, calculate feature strength by correlating the source image with the simulated linear phase band, if false (faster processing)
* these two fields will be returned equal 0.0
* @param debugLevel Debug level, turning on more text/graphic output when higher
* @return array with the following result fields (or null if feature failed to be detected):
* 0 - absolute strength of the feature (proportional to the image contrast)
* 1 - Strength reference level (RMS of the frequency samples amplitudes (in the selected band in the frequency domain)
* 2 - phaseStrength - composite value calculated when selecting the feature direction - this value is calculated even if calculateStrength is false
* 3 - angle from the tile center to the detected feature (perpendicular to the feature direction itself), zero - left, pi/2 - down (pixel coordiantes)
* 4 - distance in pixels from the tile center to the feature
* 5 - phaseAtZero Detected feature type (may be constrained by featureFilter - see that parameter for meaning of the different phase values
*/
public double [] tileLineDetect(
DoubleFHT doubleFHT,
int imageNumber, //boolean side,
double alphaThreshold,
final boolean useBinaryAlpha,
int corrFFTSize,
int corrXC,
int corrYC,
double correlationHighPassSigma,
double correlationLowPassSigma,
double phaseIntegrationWidth,
double resultHighPass, //cutting "ribbon" near zero frequency - influences length of the detected line (used only for strength calculation)
double dispertionCost,
int featureFilter,
double minRMSFrac,
double minAbs,
double maxPhaseMismatch,
boolean calculateStrength,
int debugLevel){
if (doubleFHT==null )doubleFHT= new DoubleFHT();
int length=corrFFTSize*corrFFTSize;
int numSensors=this.channel.length;
int numLayers=this.overlapImages.length/numSensors;
double [][] selection = getSelection(
0, // cross-correlation always
corrFFTSize,
corrXC,
corrYC,
0); // positive - shift second image, negative - shift first image
double[] hamming1d=doubleFHT.getHamming1d(corrFFTSize,0.0); // complete zero out on the border
double [] window=new double[length]; // TODO:cache window in doubleFHT instance?
for (int iy=0;iy0.0){
double sumAlpha=0.0;
double sumWindow=0.0;
for (int index=0;index0.0) || (correlationLowPassSigma>0.0)){
frequencyFilter=doubleFHT.createFrequencyFilter(
tileData, // just for length
correlationHighPassSigma,
correlationLowPassSigma);
}
//calculateAmplitude
doubleFHT.swapQuadrants(tileData);
doubleFHT.transform(tileData,false);
if (frequencyFilter!=null){
doubleFHT.multiplyByReal(tileData, frequencyFilter);
}
doubleFHT.debug=(debugLevel>3);
double [] dirDistStrength= doubleFHT.calcPhaseApproximation(
phaseIntegrationWidth,
tileData,
minRMSFrac,
minAbs,
maxPhaseMismatch,
dispertionCost,
featureFilter);
if (dirDistStrength==null) return null;
double angle= dirDistStrength[0];
double distance= dirDistStrength[1];
double phaseAtZero= dirDistStrength[2];
double phaseStrength=dirDistStrength[3];
if (distance<0){
angle+=Math.PI;
angle+=(2*Math.PI)*Math.round(angle/(2*Math.PI));
phaseAtZero=-phaseAtZero;
distance=-distance;
}
double [] result= {
Double.NaN, //strengths[0], // absolute strength
Double.NaN, //strengths[1], // reference level
phaseStrength, // phase strength (less computations than the first two, known anyway)
angle,
distance,
phaseAtZero
};
if (!calculateStrength) return result;
// needed data is already received, next just for debugging. Not only - also to measure the strength
int [] bandIndices= doubleFHT.getBandIndices(
false,
phaseIntegrationWidth,
angle);
double [] bandMask= doubleFHT.getRibbonMask(
false,
bandIndices,
phaseIntegrationWidth,
resultHighPass,
angle);
double [] phaseStepCorr=doubleFHT.compensatePhaseStep(
angle,
phaseAtZero,
distance,
0, //phaseTolerance, // if >0 will zero amplitude if the phase is too off - - only used with non-zero modifiedAmpPhase
null, // modifiedAmpPhase, // modifies phase and/or amplitude (if phaseTolerance>0),
bandIndices,
0); // zeroBinHalfSize); - only used with non-zero modifiedAmpPhase
double [] strengths=doubleFHT.linearFeatureStrength(
tileData,
angle,
distance,
phaseIntegrationWidth,
bandIndices,
null, // modifiedAmpPhase,
phaseStepCorr,
bandMask);
result[0]=strengths[0]; // absolute strength
result[1]=strengths[1]; // reference level
if (debugLevel>2){
System.out.println("balanced strength="+IJ.d2s(Math.sqrt(strengths[0]*strengths[0]/strengths[1]),3));
System.out.println("absolute strength="+IJ.d2s(strengths[0],3));
System.out.println("relative strength="+IJ.d2s((strengths[0]/strengths[1]),3));
System.out.println(" phase strength="+IJ.d2s(phaseStrength,3));
System.out.println(" Angle="+IJ.d2s(180*angle/Math.PI,2)+" degrees");
System.out.println(" distance="+IJ.d2s(distance,2)+" pixels");
System.out.println(" zero phase="+IJ.d2s(180*phaseAtZero/Math.PI,2)+" degrees");
}
return result;
}
public void testLineDetect(
int corrFFTSize,
int corrXC,
int corrYC,
double phaseCorrelationFraction,
double correlationHighPassSigma,
double correlationLowPassSigma,
boolean removeIslands,
double threshold, // relative to RMS value
boolean adjustDirection,
double phaseIntegrationWidth,
double resultHighPass,
double dispertionCost,
int featureFilter,
double zeroBinHalfSize,
double minRMSFrac,
double minAbs,
double maxPhaseMismatch,
double phaseTolerance, //simulation phase tolerance
int debugLevel){
int length=corrFFTSize*corrFFTSize;
DoubleFHT doubleFHT= new DoubleFHT();
int numLayers=this.overlapImages.length/2;
double [][] selection = getSelection(
1, // autocorrelation always
corrFFTSize,
corrXC,
corrYC,
0); // positive - shift second image, negative - shift first image
double[] hamming1d=(new DoubleFHT()).getHamming1d(corrFFTSize,0.0); // complete zero out on the border
double [] window=new double[length];
for (int iy=0;iy2){
(new ShowDoubleFloatArrays()).showArrays(
first,
corrFFTSize,
corrFFTSize,
"windowed-Y"+corrXC+"-y"+corrYC+"-PC"+phaseCorrelationFraction);
}
double [] frequencyFilter=null;
if ((correlationHighPassSigma>0.0) || (correlationLowPassSigma>0.0)){
frequencyFilter=doubleFHT.createFrequencyFilter(
first, // just for length
correlationHighPassSigma,
correlationLowPassSigma);
}
//calculateAmplitude
doubleFHT.swapQuadrants(first);
doubleFHT.transform(first,false);
if (frequencyFilter!=null){
doubleFHT.multiplyByReal(first, frequencyFilter);
if (debugLevel>2){
double [] filteredAmplitude=doubleFHT.calculateAmplitude(first);
(new ShowDoubleFloatArrays()).showArrays(
filteredAmplitude,
corrFFTSize,
corrFFTSize,
"filteredAmp-Y"+corrXC+"-y"+corrYC+"-PC"+phaseCorrelationFraction);
for (int i=0;i1);
double [] dirDistStrength= doubleFHT.calcPhaseApproximation(
phaseIntegrationWidth,
first,
minRMSFrac,
minAbs,
maxPhaseMismatch,
dispertionCost,
featureFilter);
if (dirDistStrength==null) return;
double angle= dirDistStrength[0];
double distance= dirDistStrength[1];
double phaseAtZero= dirDistStrength[2];
double phaseStrength=dirDistStrength[3];
// needed data is already received, next just for debugging
int [] bandIndices= doubleFHT.getBandIndices(
false,
phaseIntegrationWidth,
angle);
double [] bandMask= doubleFHT.getRibbonMask(
false,
bandIndices,
phaseIntegrationWidth,
resultHighPass,
angle);
// double [][] modifiedAmpPhase=doubleFHT.calcIndAmHPhase (first, bandIndices);
double [] phaseStepCorr=doubleFHT.compensatePhaseStep(
angle,
phaseAtZero,
distance,
phaseTolerance, // if >0 will zero amplitude if the phase is too off
null, // modifiedAmpPhase, // modifies phase and/or amplitude (if phaseTolerance>0),
bandIndices,
zeroBinHalfSize);
double [] strengths=doubleFHT.linearFeatureStrength(
first,
angle,
distance,
phaseIntegrationWidth,
bandIndices,
null, // modifiedAmpPhase,
phaseStepCorr,
bandMask);
if (debugLevel>0){
System.out.println("balanced strength="+IJ.d2s(Math.sqrt(strengths[0]*strengths[0]/strengths[1]),3));
System.out.println("absolute strength="+IJ.d2s(strengths[0],3));
System.out.println("relative strength="+IJ.d2s((strengths[0]/strengths[1]),3));
System.out.println(" phase strength="+IJ.d2s(phaseStrength,3));
System.out.println(" Angle="+IJ.d2s(180*angle/Math.PI,2)+" degrees");
System.out.println(" distance="+IJ.d2s(distance,2)+" pixels");
System.out.println(" zero phase="+IJ.d2s(180*phaseAtZero/Math.PI,2)+" degrees");
}
if (debugLevel>0){
double [] reconstructed=doubleFHT.reconstructLinearFeature(
ones,
angle,
distance,
phaseIntegrationWidth,
bandIndices,
null, //modPhase,
phaseStepCorr,
bandMask);
//debugWindowed
double [][] reconstructionComparison={reconstructed, debugWindowed};
String [] titles={"reconstructed","windowed original"};
(new ShowDoubleFloatArrays()).showArrays(
reconstructionComparison,
corrFFTSize,
corrFFTSize,
true,
"reconstructed-x"+corrXC+"-y"+corrYC,
titles);
}
}
public void testLineDetectDebug(
int corrFFTSize,
int corrXC,
int corrYC,
double phaseCorrelationFraction,
double correlationHighPassSigma,
double correlationLowPassSigma,
boolean removeIslands,
double threshold, // relative to RMS value
boolean adjustDirection,
double phaseIntegrationWidth,
double resultHighPass,
double dispertionCost,
int featureFilter,
double zeroBinHalfSize,
double minRMSFrac,
double minAbs,
double maxPhaseMismatch,
double phaseTolerance, //simulation phase tolerance
int debugLevel){
int length=corrFFTSize*corrFFTSize;
DoubleFHT doubleFHT= new DoubleFHT();
// double [][] selection = new double[this.overlapImages.length][length];
int numLayers=this.overlapImages.length/2;
double [][] selection = getSelection(
1, // autocorrelation always
corrFFTSize,
corrXC,
corrYC,
0); // positive - shift second image, negative - shift first image
// correlation stuff
// double[] hamming1d=(new DoubleFHT()).getHamming1d(corrFFTSize);
double[] hamming1d=(new DoubleFHT()).getHamming1d(corrFFTSize,0.0); // complete zero out on the border
double [] window=new double[length];
for (int iy=0;iy2){
(new ShowDoubleFloatArrays()).showArrays(
first,
corrFFTSize,
corrFFTSize,
"windowed-Y"+corrXC+"-y"+corrYC+"-PC"+phaseCorrelationFraction);
}
double [] frequencyFilter=null;
if ((correlationHighPassSigma>0.0) || (correlationLowPassSigma>0.0)){
frequencyFilter=doubleFHT.createFrequencyFilter(
first, // just for length
correlationHighPassSigma,
correlationLowPassSigma);
}
//calculateAmplitude
doubleFHT.swapQuadrants(first);
doubleFHT.transform(first,false);
if (debugLevel>2){
double [] unfilteredAmplitude=doubleFHT.calculateAmplitude(first);
for (int i=0;i2){
double [] filteredAmplitude=doubleFHT.calculateAmplitude(first);
(new ShowDoubleFloatArrays()).showArrays(
filteredAmplitude,
corrFFTSize,
corrFFTSize,
"filteredAmp-Y"+corrXC+"-y"+corrYC+"-PC"+phaseCorrelationFraction);
for (int i=0;i1);
double [] dirDistStrength= doubleFHT.calcPhaseApproximation(
phaseIntegrationWidth,
first,
minRMSFrac,
minAbs,
maxPhaseMismatch,
dispertionCost,
featureFilter);
if (dirDistStrength==null) return;
double angle= dirDistStrength[0];
double distance= dirDistStrength[1];
double phaseAtZero= dirDistStrength[2];
// double phaseStrength=dirDistStrength[3];
// needed data is already received, next just for debugging
int [] bandIndices= doubleFHT.getBandIndices(
false,
phaseIntegrationWidth,
angle);
double [] bandMask= doubleFHT.getRibbonMask(
false,
bandIndices,
phaseIntegrationWidth,
resultHighPass,
angle);
double [][] modifiedAmpPhase=doubleFHT.calcIndAmHPhase (first, bandIndices);
double [] phaseStepCorr=doubleFHT.compensatePhaseStep(
angle,
phaseAtZero,
distance,
phaseTolerance, // if >0 will zero amplitude if the phase is too off
modifiedAmpPhase, // modifies phase and/or amplitude (if phaseTolerance>0),
bandIndices,
zeroBinHalfSize);
if (debugLevel>0){
double [][] dbgMask=new double [2][corrFFTSize*corrFFTSize/2];
for (int i=0;i0){
System.out.println(" ====== Linear feature strength: absolute="+strengths[0]+" relative="+(strengths[0]/strengths[1]));
}
}
public double [] correlateRectangular(
boolean autoCorrelation,
int corrFFTSize,
int overlapFraction,
int corrXC,
int corrYC,
int maxDisparity,
double phaseCorrelationFraction,
double corrCbWeight,
double corrCrWeight,
double correlationHighPassSigma,
double correlationLowPassSigma,
double noiseNormalizationSignaY,
double noiseNormalizationSignaCbCr,
double contrastThreshold,
boolean useBinaryAlpha,
boolean enableNegativeDisparity,
int threadsMax,
int debugLevel){
int numTiles=1+(maxDisparity/(corrFFTSize/overlapFraction));
int length=corrFFTSize*corrFFTSize;
int corrWidth=corrFFTSize + corrFFTSize*(numTiles-1)/overlapFraction;
int outLength=corrFFTSize*corrWidth;
double [] fullCorr=new double [outLength];
for (int i=0;i0){
System.out.println("Selection ("+(useBinaryAlpha?"binary":"analog")+" alpha ) xc="+corrXC+" yc="+corrYC+" size="+corrFFTSize+" maxDisparity="+maxDisparity+
" : first "+IJ.d2s(100*stats[0],3)+"%, second "+IJ.d2s(100*stats[1],3)+"%, overlap "+IJ.d2s(100*stats[2],3)+"%");
}
double [][] window={selection[0].clone(),selection[numLayers].clone()};
for (int iy=0;iy1)){
debugCorr[1]=first.clone();
debugCorr[2]=second.clone();
}
// TODO: use common (per-thread) DoubleFHT instance to use caching of sin,cos, filter tables
corr[i+1]=(new DoubleFHT()).correlate (
first,
second,
correlationHighPassSigma,
correlationLowPassSigma,
phaseCorrelationFraction);
double sigma=(i==0)?noiseNormalizationSignaY:noiseNormalizationSignaCbCr;
if (sigma>0){
double rms=hFNoiseNormalize(
corr[i+1], // double [] data,
sigma, // double filterSigma,
true); // boolean centerOnly);
if (this.debugLevel>1){
System.out.println("Correlation RMS for component "+((i==0)?"Y":((i==1)?"Cb":"Cr"))+ " was "+rms);
}
}
// for (int j=0;j2){
String [] titles={"corr","first","second"};
debugCorr[0]=corr[0];
(new ShowDoubleFloatArrays()).showArrays(
debugCorr,
corrFFTSize,
corrFFTSize,
true,
"corr"+nTile+"_selection-x"+corrXC+"-y"+corrYC+"-MAXDSP"+maxDisparity+"-PC"+phaseCorrelationFraction,
titles);
}
int srcIndex=0;
for (int y=0;y>=1) if ((chm & 1)!=0) numLayers++;
if (externalData!=null) numLayers++;
int [] img= {first,second};
int [] channels=new int [numLayers];
int c=0;
int l=0;
for (int chm=channelMask; chm!=0; chm>>=1) {
if ((chm & 1)!=0) channels[l++]=c;
c++;
}
if (externalData!=null) channels[l]=3;
if (debugLevel>1){
double d=Math.sqrt(dxA*dxA+dyA*dyA);
System.out.println(" getShiftedSlices(): xc1="+xc1+" yc1="+yc1+" size="+size+" first="+first+" second="+second);
System.out.println(" getShiftedSlices(): d="+IJ.d2s(d,2)+" dX="+IJ.d2s(dxA,2)+" dY="+IJ.d2s(dyA,2));
System.out.println(" getShiftedSlices(): shiftImage[0]="+shiftImage[0]+" shiftImage[1]="+shiftImage[1]);
System.out.println(" getShiftedSlices() first image: dX="+IJ.d2s(dx[0],2)+" dY="+IJ.d2s(dy[0],2));
System.out.println(" getShiftedSlices() second image: dX="+IJ.d2s(dx[1],2)+" dY="+IJ.d2s(dy[1],2));
System.out.println(" getShiftedSlices() second image: iDx="+iDx+" iDyY="+iDy);
System.out.println(" getShiftedSlices() second image: xc="+(xc1+iDx)+" yc="+(yc1+iDy));
}
int doubleSize=2*size;
int doubleLength=doubleSize*doubleSize;
double [] slice= new double[doubleLength];
double [][] result=new double [2*numLayers][];
double [] imgSlice;
if (window==null) {
window = new double [doubleLength];
double[] window1d=doubleFHT.getHamming1d(size); // Hamming
int index=0;
int halfSize=size/2;
int size32=3*size/2;
for (int iy=0;iysize32)?window1d[iy-size]:1.0);
for (int ix=0;ixsize32)?window1d[ix-size]:1.0);
window[index++]=wx*wy;
}
}
}
if (debugLevel>1){
// double d=Math.sqrt(dxA*dxA+dyA*dyA);
System.out.print("\ngetShiftedSlices(): channelMask="+channelMask+" channels=");
for (int i=0;i2){
(new ShowDoubleFloatArrays()).showArrays(
window,
doubleSize,
doubleSize,
"window");
}
}
for (int iImg=0;iImg=0) && (chn<3)){ // starts from 0 (Y)
imgSlice=this.overlapImages[numImgChannels*nImg+chn+1]; // skip alpha
if (debugLevel>2) System.out.println("Using slice "+(numImgChannels*nImg+chn+1));
} else {
imgSlice=externalData[nImg];
if (debugLevel>2) System.out.println("Using externalData["+nImg+"]");
}
if (debugLevel>2) System.out.println("xc["+iImg+"]="+xc[iImg]+" yc["+iImg+"]="+yc[iImg]+" shiftImage["+iImg+"]="+shiftImage[iImg]);
slice= getSelection(
imgSlice, // source image/channel slice
slice,
doubleSize, //int width,
doubleSize, //int height,
xc[iImg] , //xc,
yc[iImg]); //yc);
if (debugLevel>2) {
(new ShowDoubleFloatArrays()).showArrays(
slice,
doubleSize,
doubleSize,
"slice");
}
// normalize and save DC
double dc=0.0;
if (shiftImage[iImg]){
dc=normalizeAndWindowGetDC (slice, window); //windowInterpolation
// doubleFHT.shift(slice, dx[iImg], dy[iImg]);
doubleFHT.shift(slice, -dx[iImg], -dy[iImg]);
}
int oSliceNumber=img.length*cN+iImg;
if (doubleSizeOutput) {
if (debugLevel>2) System.out.println("doubleLength="+doubleLength+"doubleSize="+doubleSize);
result[oSliceNumber]=new double [doubleLength];
int oIndex=0;
for (int iY=0;iY=result[oSliceNumber].length) || (oIndex>=slice.length) || (oIndex>=window.length)){
System.out.println("iX="+iX+" iY="+iY+" oIndex="+oIndex);
}
// result[oSliceNumber][oIndex++]=imgSlice[oIndex]/window[oIndex]+dc;// no NaN with Hamming
result[oSliceNumber][oIndex]=slice[oIndex]/window[oIndex]+dc;// no NaN with Hamming
oIndex++;
}
}
}else{
result[oSliceNumber]=new double [size*size];
int iIndex=size*size+size/2; // top left index of the inner square size*size
int oIndex=0;
for (int iY=0;iY2){
String [] channelNames={"Y","Cb","Cr","Lin"};
String [] titles=new String[2*numLayers];
String title="SH";
for (int i=0;i=nImg) nSImg++; // here - second image number, not index!
imp_stack.setProperty("firstImage_"+nPair, ""+nImg);
imp_stack.setProperty("secondImage_"+nPair, ""+nSImg);
}
if (this.disparityScales!=null) {
imp_stack.setProperty("numberOfImages", this.disparityScales.length+"");
for (int i=0;ipairs[nPair][0])?(pairs[nPair][1]-1):pairs[nPair][1];
}
this.corrFFTSize= Integer.parseInt ((String) imp.getProperty("corrFFTSize"));
this.overlapStep= Integer.parseInt ((String) imp.getProperty("overlapStep"));
this.disparityPerEntry= Double.parseDouble ((String) imp.getProperty("disparityPerEntry"));
this.tilesY= Integer.parseInt ((String) imp.getProperty("tilesY"));
this.tilesX= Integer.parseInt ((String) imp.getProperty("tilesX"));
if (imp.getProperty("numberOfImages")==null){ // later - add to required properties, now just to use older files
int numImages=-1;
for (int i=0;inumImages) numImages=pairs[i][0];
if (pairs[i][1]>numImages) numImages=pairs[i][1];
}
numImages++;
double [][] disparityPair={{-0.5,0.0},{0.5,0.0}};
double [][] disparityTriplet={{0.0,Math.sqrt(3)/3}, {-0.5,-Math.sqrt(3)/6},{0.5,-Math.sqrt(3)/6}};
if ((numImages==2)|| (numImages==3)){
this.disparityScales=new double [numImages][2];
for (int i=0;i2 ) && (tile== ( (tiles/tilesX/2)*tilesX +(tilesX/2)));
if (debugThis){
System.out.println("tile="+tile+" tiles="+tiles+" tilesX="+tilesX);
}
for (int nPair=0;nPair2){
(new ShowDoubleFloatArrays()).showArrays(
this.centerPixels,
this.impDisparity.getWidth(),
this.centerPixels.length/this.impDisparity.getWidth(),
"center-corr");
}
}
/**
* Calculating combined correlation arrays for the actual cameras from the center camera
* Will do nothing if the data is already calcualted
* @param threadsMax maximal number of therads to use
* @param showProgress show ImageJ progress/status
* @param debugLevel debug level
*/
public void moveDisparityFromCenter(
final int threadsMax,
final boolean showProgress,
final int debugLevel){
if (this.syntheticPixels!=null) return; // already done
if (this.centerPixels==null) {
String msg="Centered disparity data is not defined";
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
final int numImages= this.disparityScales.length;
this.syntheticPixels=new float [numImages] [this.centerPixels.length]; // all should be the same
// TODO - convert to multithread?
final int tilesPerImage=this.tilesY*this.tilesX;
final int tiles=tilesPerImage*numImages;
final int tilesX=this.tilesX;
final double [][] disparityScales=this.disparityScales;
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger tileIndexAtomic = new AtomicInteger(0);
final float [] centerPixels=this.centerPixels;
final float [][] syntheticPixels=this.syntheticPixels;
if (showProgress) IJ.showProgress(0.0);
if (showProgress) IJ.showStatus("Re-centering correlation data...");
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
@Override
public void run() {
for (int tile=tileIndexAtomic.getAndIncrement(); tile2 ) && (nImage==0)&& (tile== ( (tiles/tilesX/2)*tilesX +(tilesX/2)));
if (debugThis){
System.out.println("nImage="+nImage+"tile="+tile+" tiles="+tiles+" tilesX="+tilesX);
}
double [] dXY={-disparityScales[nImage][0],-disparityScales[nImage][1]};
if (debugThis){
System.out.println("\n\nnImg="+nImage+" dXY[0]="+dXY[0]+" dXY[1]="+dXY[1]);
}
double [] tileData=getRecenteredDisparity(
tileX, //int tileX,
tileY, //int tileY,
centerPixels,
dXY,
debugThis?(debugLevel+2):debugLevel);
if (tileData==null) System.out.println("combineDisparityToCenter(): tileData==null, tile="+tile);
setTileDisparity(
tileData,
tileX,
tileY,
syntheticPixels[nImage]);
if (showProgress){
final int fTile=tile;
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
// Here, we can safely update the GUI
// because we'll be called from the
// event dispatch thread
IJ.showProgress(fTile,tiles);
}
});
}
}
}
};
}
startAndJoin(threads);
if (showProgress) IJ.showStatus("");
if (showProgress) IJ.showProgress(1.0);
}
private void initDoubleTileWindow(){
int size=this.overlapStep*2;
this.doubleTileWindow = new double [size*size];
double [] window1d=(new DoubleFHT()).getHamming1d(size,0); // pure cosine
for (int i=0;i3) System.out.println("getRecenteredDisparity("+tileX+","+tileY+",pixels,...)");
int disparityPoints=getDisparityPoints();
double disparityPerEntry=getDisparityPerEntry();
double lDisp=Math.sqrt(relDxy[0]*relDxy[0]+relDxy[1]*relDxy[1])*(disparityPoints-1)*disparityPerEntry; // in Pixels
int maxTileDis=(int) Math.floor(lDisp/this.overlapStep)+2; // 1 is enough, maybe some rounding?
double [][][] tileCache=new double [2*maxTileDis+1][2*maxTileDis+1][];
for (int i=0;i3) System.out.print(" maxTileDis="+maxTileDis);
for (int i=0;i3) System.out.print("getRecenteredDisparity() i="+i);
double dTX=kTXY[0]*i; // tile period is 0.5 in this units, all tiles within -1.0...+1.0 may contribute - normally 4
double dTY=kTXY[1]*i;
double d=0.0;
double weight=0.0;
int idTX0=(int) Math.floor(dTX);
int idTY0=(int) Math.floor(dTY);
for (int iTY=0;iTY<2;iTY++) {
int iTYCache=idTY0+iTY+maxTileDis;
int pY= (int) Math.round ((dTY-(idTY0+iTY)+1.0)*this.overlapStep); // 0... 2*this.overlapStep
int srcTY=tileY+idTY0+iTY;
if ((srcTY>=0) && (srcTY=0) && (pY=0) && (srcTX=0) && (pX3) || ((tileX+idTX0+iTX)<0) || ((tileY+idTY0+iTY)<0)){
System.out.println ("\n>>>getting tile to cache: iTXCache="+iTXCache+" iTYCache="+iTYCache+" tileX="+tileX+" tileY="+tileY+" idTX0="+idTX0+" idTY0="+idTY0+ " iTX="+iTX+" iTY="+iTY);
}
tileCache[iTYCache][iTXCache]=getTileDisparity( // negative?
srcTX,
srcTY,
tiledCorrelation);
}
double w=this.doubleTileWindow[pY*windowSize+pX];
weight+=w;
d+=w*tileCache[iTYCache][iTXCache][i];
if (debugLevel>3) {
System.out.print(" pX="+pX+" pY="+pY+" srcTX="+srcTX+" srcTY="+srcTY+" d="+tileCache[iTYCache][iTXCache][i]+" w="+w);
}
}
}
}
result[i]=(weight>0)?(d/weight):0.0;
if (debugLevel>3) System.out.println(" weight="+weight+" result["+i+"]="+result[i]);
}
if (debugLevel>3) {
for (int i=0;i=data.length) || (i>=disparityArray.length)){
System.out.println("setTileDisparity(disparityArray,"+itx+","+ity+", data): index="+index+
" data.length="+data.length+" i="+i+" disparityArray.length="+disparityArray.length+" disparityPoints="+disparityPoints+" width="+width);
}
data[index++]=(float) disparityArray[i]; //oob 20375037
}
}
public double [] getCorrelationArray(
double x,
double y,
double fatZero,
int threadsMax,
boolean showProgress,
int debugLevel ){
combineDisparityToCenter(
fatZero,
threadsMax,
showProgress,
debugLevel);
return getCorrelationArray(
this.centerPixels,
x,
y,
threadsMax,
showProgress,
debugLevel );
}
public double [] getCorrelationArray(
int imageNumber,
double x,
double y,
int threadsMax,
boolean showProgress,
int debugLevel ){
moveDisparityFromCenter(
threadsMax,
showProgress,
debugLevel);
float [] data=this.syntheticPixels[imageNumber];
return getCorrelationArray(
data,
x,
y,
threadsMax,
showProgress,
debugLevel );
}
public double [] getCorrelationArray(
int firstImage,
int secondImageIndex,
// double fatZero,
double x,
double y,
int threadsMax,
boolean showProgress,
int debugLevel ){
int nPair;
nPair= imagesToNPair(firstImage, secondImageIndex);
if (nPair<0) {
System.out.println("getCorrelationArray(): No data for pair of images "+firstImage+" and second image index "+secondImageIndex+
" (second image = "+((secondImageIndex>=firstImage)?(secondImageIndex+1):(secondImageIndex))+")");
return null;
}
float [] data=(nPair>=0)?this.pixels[nPair]:this.centerPixels;
return getCorrelationArray(
data,
x,
y,
threadsMax,
showProgress,
debugLevel );
}
private double [] getCorrelationArray(
float [] data,
double x,
double y,
int threadsMax,
boolean showProgress,
int debugLevel ){
double tx=x/this.overlapStep;
double ty=y/this.overlapStep;
int itx0=(int) Math.floor(tx);
int ity0=(int) Math.floor(ty);
tx-=itx0;
ty-=ity0;
int itx1=itx0+1;
int ity1=ity0+1;
if ((itx1<0) || (itx0>=this.tilesX) || (ity1<0) || (ity0>=this.tilesY)) {
System.out.println("getCorrelationArray(): this.overlapStep="+this.overlapStep+" this.tilesX="+this.tilesX+"this.tilesY="+this.tilesY);
System.out.println("getCorrelationArray(): itx0="+itx0+" ity0="+ity0+" out of bounds. x="+x+" y="+y);
return null;
}
if (itx0<0){itx0=itx1;tx=0; }
if (itx1>=this.tilesX){itx1=itx0;tx=0;}
if (ity0<0){ity0=ity1;ty=0; }
if (ity1>=this.tilesY){ity1=ity0;ty=0;}
int width=this.impDisparity.getWidth();
int disparityPoints=width/this.tilesX;
int index00=ity0*width+itx0*disparityPoints;
int index10=ity0*width+itx1*disparityPoints;
int index01=ity1*width+itx0*disparityPoints;
int index11=ity1*width+itx1*disparityPoints;
double [] result=new double [disparityPoints];
for (int i=0;iresult.length) ||
((index00+i)>=data.length) || ((index10+i)>=data.length) || ((index01+i)>=data.length) || ((index11+i)>=data.length)){
System.out.println("width="+width+" index00="+index00+" index10="+index10+" index01="+index01+" index11="+index11);
System.out.println("itx0="+itx0+" itx1="+itx1+" ity0="+ity0+" ity1="+ity1+" this.tilesX="+this.tilesX+" this.tilesY="+this.tilesY);
}
result[i]= // bi-linear interpolation //TODO:java.lang.ArrayIndexOutOfBoundsException: 24760512at PixelMapping$InterSensor$DisparityTiles.getCorrelationArray(PixelMapping.java:4619)
(data[index00+i]*(1.0-tx)+data[index10+i]*tx)*(1.0-ty)+
(data[index01+i]*(1.0-tx)+data[index11+i]*tx)*ty;
}
return result;
}
public double [] getCorrelationArray(
int imageNumber,
double fatZero,
double x,
double y,
int threadsMax,
boolean showProgress,
int debugLevel ){
int numSecondImages=0;
for (int i=0;i0) d*=(partial[nI][nP]+fatZero);
else d=0;
}
result[nP]=Math.pow(d, pwr)-fatZero;
}
return result;
}
private int imagesToNPair(int firstImage, int secondImageIndex){
for (int i=0;i=0)?woi.x:0;
this.zMapWOI.y=(woi.y>=0)?woi.y:0;
this.zMapWOI.width= (((woi.x+woi.width) <=this.tilesX)?(woi.x+woi.width): this.tilesX)-this.zMapWOI.x;
this.zMapWOI.height=(((woi.y+woi.height)<=this.tilesY)?(woi.y+woi.height):this.tilesY)-this.zMapWOI.y;
final Rectangle zMapWOI=this.zMapWOI;
if ((this.zMap==null) || (this.zMap[nImg]==null)) initZMap(nImg);
final ZTile [][] thisZMap=this.zMap[nImg];
final int tilesX=this.zMapWOI.width;
final int tiles=this.zMapWOI.width*this.zMapWOI.height;
if (debugLevel>2) System.out.println("fillForegroundGaps() woi.x="+woi.x+" woi.y="+woi.y+
" woi.width="+woi.width+" woi.height="+woi.height);
if (debugLevel>2) System.out.println("fillForegroundGaps() zMapWOI.x="+zMapWOI.x+" zMapWOI.y="+zMapWOI.y+
" zMapWOI.width="+zMapWOI.width+" zMapWOI.height="+zMapWOI.height);
if (debugLevel>1) System.out.println("fillForegroundGapsStep() foregroundThreshold="+foregroundThreshold+" maxDifference="+maxDifference+" tiles="+tiles);
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger tileIndexAtomic = new AtomicInteger(0);
if (showProgress) IJ.showProgress(0.0);
if (showProgress) IJ.showStatus("Setting up zMap for image "+nImg+" ("+(nImg+1)+" of "+this.disparityScales.length+") ...");
final AtomicInteger numberUpdated = new AtomicInteger(0);
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
@Override
public void run() {
for (int tile=tileIndexAtomic.getAndIncrement(); tile4) System.out.println("0:fillForegroundTile("+tileX+","+tileY+",zMap,"+minNeib+","+maxDifference+","+foregroundThreshold+","+debugLevel+"), zTile.foregroundIndex="+zTile.foregroundIndex);
if (zTile.foregroundIndex==0) return false; // already front
if (debugLevel>3) System.out.println("fillForegroundTile("+tileX+","+tileY+",zMap,"+minNeib+","+maxDifference+","+foregroundThreshold+","+debugLevel+"), zTile.foregroundIndex="+zTile.foregroundIndex);
int [][] dirs8={{1,0},{1,1},{0,1},{-1,1},{-1,0},{-1,-1},{0,-1},{1,-1}};
for (int plane=0; plane =foregroundThreshold){
// See if there are enough neighbors where current foreground is close enough to this one
int numMatched=0;
for (int iDir=0;iDir=zMap.length){
// System.out.print("fillForegroundTile("+tileX+","+tileX+",...) zMap.length="+zMap.length+" zMap[0].length="+zMap[0].length);
// }
if ((tileY1>=0) && (tileY1=0) && (tileX1=minNeib){
zTile.setForegroundPlane(plane);
return true;
}
}
return false;
}
public void foregroundByOcclusion(
double [][][] imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
int imageFullWidth,
double bgFraction, // process maximal correlation or the first farther plane with correlation intensity >= this part of maximal
double blurVarianceSigma,
int refineTilePeriod,
double refinePhaseCoeff,
double refineHighPassSigma,
double refineLowPassSigma,
double refineCorrMaxDistance,
double refineCorrThreshold,
int refineSubPixel,
double zMapMinForeground,
int zMapVarMask,
double [] zMapVarThresholds,
double [] zMapVarWeights,
int normalizeCorrelation,
int zMapCorrMask,
double [] zMapCorrThresholds,
double [] zMapCorrWeights,
int debugRow,
int debugColumn,
int threadsMax,
boolean showProgress,
int debugLevel){
// process all image pairs
for (int nImg=0;nImg= this part of maximal
blurVarianceSigma,
refineTilePeriod,
refinePhaseCoeff,
refineHighPassSigma,
refineLowPassSigma,
refineCorrMaxDistance,
refineCorrThreshold,
refineSubPixel,
zMapMinForeground,
zMapVarMask,
zMapVarThresholds,
zMapVarWeights,
normalizeCorrelation,
zMapCorrMask,
zMapCorrThresholds,
zMapCorrWeights,
debugRow,
debugColumn,
threadsMax,
showProgress,
debugLevel);
}
}
public void foregroundByOcclusion( // remove unneeded parameters later
final int nImg,
final int[] sImg, // list of second images
final double [][][] imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
final int imageFullWidth,
final double bgFraction, // process maximal correlation or the first farther plane with correlation intensity >= this part of maximal
final double blurVarianceSigma,
final int refineTilePeriod,
final double refinePhaseCoeff,
final double refineHighPassSigma,
final double refineLowPassSigma,
final double refineCorrMaxDistance,
final double refineCorrThreshold,
final int refineSubPixel,
final double zMapMinForeground,
final int zMapVarMask,
final double [] zMapVarThresholds,
final double [] zMapVarWeights,
final int normalizeCorrelation,
final int zMapCorrMask,
final double [] zMapCorrThresholds,
final double [] zMapCorrWeights,
final int debugRow,
final int debugColumn,
final int threadsMax,
final boolean showProgress,
final int debugLevel){
final int debugThreshold=2;
final Rectangle zMapWOI=this.zMapWOI;
final ZTile [][] thisZMap=this.zMap[nImg];
final int tilesX=this.zMapWOI.width;
final int tiles=this.zMapWOI.width*this.zMapWOI.height;
if (debugLevel>2) System.out.println("setupZMap() zMapWOI.x="+zMapWOI.x+" zMapWOI.y="+zMapWOI.y+
" zMapWOI.width="+zMapWOI.width+" zMapWOI.height="+zMapWOI.height);
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger tileIndexAtomic = new AtomicInteger(0);
if (showProgress) IJ.showProgress(0.0);
if (showProgress) IJ.showStatus("Filtering zMap foreground for image "+nImg+" ("+(nImg+1)+" of "+this.disparityScales.length+") ...");
final int debugTile=debugRow*tilesX+debugColumn;
//TODO: define here
final double [] window=new double [4*this.overlapStep*this.overlapStep];
window[0]=Double.NaN;
final double [] refineWindow=new double [refineTilePeriod*refineTilePeriod*4];
refineWindow[0]=Double.NaN;
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
@Override
public void run() {
DoubleFHT doubleFHT = new DoubleFHT();
DoubleFHT refineFHT= new DoubleFHT();
for (int tile=tileIndexAtomic.getAndIncrement(); tile=debugThreshold) && (tile==debugTile);
int thisDebugLevel=debugLevel+(debugThis?2:0);
// while (!
foregroundByOcclusionTile (
zMapWOI.x+tile%tilesX, //int tileX,
zMapWOI.y+tile/tilesX, //int tileY,
thisZMap,
nImg,
sImg, // list of second images
bgFraction, // process maximal correlation or the first farther plane with correlation intensity >= this part of maximal
imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
imageFullWidth,
blurVarianceSigma,
refineTilePeriod,
refinePhaseCoeff,
refineHighPassSigma,
refineLowPassSigma,
refineCorrMaxDistance,
refineCorrThreshold,
refineSubPixel,
zMapMinForeground,
zMapVarMask,
zMapVarThresholds,
zMapVarWeights,
normalizeCorrelation,
zMapCorrMask,
zMapCorrThresholds,
zMapCorrWeights,
window,
doubleFHT,
refineWindow,
refineFHT,
threadsMax,
showProgress,
thisDebugLevel); //);
if (showProgress){
final int finalTile=tile;
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
IJ.showProgress(finalTile,tiles);
}
});
}
}
}
};
}
startAndJoin(threads);
IJ.showProgress(1.0);
}
private void foregroundByOcclusionTile ( // false if nothing left in the current foreground, may repeat
int tileX,
int tileY,
ZTile [][] thisZMap,
int nImg,
int [] sImgSet, // list of second images
double bgFraction, // process maximal correlation or the first farther plane with correlation intensity >= this part of maximal
double [][][] imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
int imageFullWidth,
double blurVarianceSigma,
int refineTilePeriod,
double refinePhaseCoeff,
double refineHighPassSigma,
double refineLowPassSigma,
double refineCorrMaxDistance,
double refineCorrThreshold,
int refineSubPixel,
double zMapMinForeground,
int zMapVarMask,
double [] zMapVarThresholds,
double [] zMapVarWeights,
int normalizeCorrelation,
int zMapCorrMask,
double [] zMapCorrThresholds,
double [] zMapCorrWeights,
double [] window,
DoubleFHT doubleFHT,
double [] refineWindow,
DoubleFHT refineFHT,
int threadsMax,
boolean showProgress,
int debugLevel){
double [] normVarWeights= normalizeWeights(zMapVarMask,zMapVarWeights);
double [] normCorrWeights=normalizeWeights(zMapCorrMask,zMapCorrWeights);
/// final double [][] disparityScales=this.disparityScales;
// this filter will iterate through second images and AND filter for foreground pixels only
ZTile zTile=thisZMap[tileY][tileX];
if (zTile==null) return; // true; // tile is not defined, nothing left to do
// find the background plane to process
int bgPlane=0;
if (zTile.maximums.length>0){
for (int i=0;izTile.maximums[bgPlane][1]) bgPlane=i;
if (bgFraction>0){
int bgPlane1=bgPlane+1;
for (int i=bgPlane+1;izTile.maximums[bgPlane1][1]) bgPlane1=i;
if ((bgPlane1=bgFraction*zTile.maximums[bgPlane][1])) bgPlane=bgPlane1;
} else if (bgFraction<0){ // compare with foreground
boolean hasFG=false;
for (int i=0;i=(-bgFraction)*zTile.maximums[bgPlane][1]){
hasFG=true;
break;
}
if (!hasFG){
bgPlane++;
if (bgPlanezTile.maximums[bgPlane][1]) bgPlane=i;
}
}
} else { // use gap-fileld FG
//zTile.foregroundIndex
if (zTile.foregroundIndex>=bgPlane){
bgPlane++;
if (bgPlanezTile.maximums[bgPlane][1]) bgPlane=i;
}
}
}
}
// now bgPlane - background plane (if >=maximums.length - use infinity)
// double bgFraction, // process maximal correlation or the first farther plane with correlation intensity >= this part of maximal
// is there any forground?
/*
if (!zTile.advanceForeground()) return true; // got to the end
if (zTile.enabledPixels[zTile.foregroundIndex]==null) {
zTile.initEnabledForegroundPixelsPlane();
}
*/
int tileOverlap=zTile.getOverlap();
int paddedSize=zTile.getPaddedSize();
int paddedLength=zTile.getPaddedLength();
// double disparity=zTile.maximums[zTile.foregroundIndex][0];
double disparity=(bgPlane>=zTile.maximums.length)?0.0:zTile.maximums[bgPlane][0];
int size=this.overlapStep*2;
int margin=size/4-tileOverlap;
int [] dirs1={1,size+1,size,size-1,-1,-size-1,-size,-size+1,0};
double [][][][] slices=new double [sImgSet.length][][][];
double [][][][] variance=new double [sImgSet.length][][][]; // [sIndex][pair 0/1][chn][pixel
int length=0; // initialize it to a correct value right here?
// int numOther=sImgSet.length;
int [][] otherPairs=new int [(sImgSet.length*(sImgSet.length-1))/2][2];
// with 3 images there is only one other pair, but try to maki it work later with more images
{
int index=0;
for (int i=0;i4){
(new ShowDoubleFloatArrays()).showArrays(
borderMask,
size,
size,
"borderMask");
}
// zTile.initAux(2);
// zTile.initAux(7);
zTile.initAux(11);
if (debugLevel>3){ // +2 for selected tile
System.out.println("otherPairs.length="+otherPairs.length+ "");
for (int i=0;i4)&& (n==0) && (chn==0) && (numPix!=9) ){
System.out.println (" x="+(i%size)+" y="+(i/size)+" i="+i+" numPix="+numPix+" ="+borderMask[i]+" v2="+v2);
}
// v2+=(S2-(S1*S1)/numPix)/numPix; //dirs1.length;
variance[sIndex][n][chn][i]=Math.sqrt(v2); //java.lang.NullPointerException
}
if (blurVarianceSigma>0.0){
if (debugLevel>3) // +2 for selected tile
System.out.println("Bluring varaince["+sIndex+"]["+n+"]["+chn+"] with sigma="+blurVarianceSigma);
(new DoubleGaussianBlur()).blurDouble(
variance[sIndex][n][chn],
size,
size,
blurVarianceSigma,
blurVarianceSigma,
0.01);
}
}
}
if (debugLevel>3){ // +2 for selected tile
int numActive=0;
for (int i=0;i4) {
System.out.print("chn="+chn+" nPair="+nPair+" i="+i+" iPix="+iPix+ " ix="+(i%paddedSize)+" iy="+(i/paddedSize)+
" ipx="+(iPix%size)+" ipy="+(iPix/size)+
" v["+sIndex0+"][0]["+chn+"]["+iPix+"]="+IJ.d2s(variance[sIndex0][0][chn][iPix],5)+
" v["+sIndex1+"][0]["+chn+"]["+iPix+"]="+IJ.d2s(variance[sIndex1][0][chn][iPix],5)+
" v["+sIndex0+"][1]["+chn+"]["+iPix+"]="+IJ.d2s(variance[sIndex0][1][chn][iPix],5)+
" v["+sIndex1+"][1]["+chn+"]["+iPix+"]="+IJ.d2s(variance[sIndex1][1][chn][iPix],5)+
" pairVar="+IJ.d2s(pairVar,5));
}
double pairDiff=Math.abs(slices[sIndex0][chn+1][1][iPix]-slices[sIndex1][chn+1][1][iPix])/pairVar/zMapVarThresholds[chn];
double thisDiff=Math.abs(0.5*(
slices[sIndex0][chn+1][0][iPix]+
slices[sIndex1][chn+1][0][iPix]-
slices[sIndex0][chn+1][1][iPix]-
slices[sIndex1][chn+1][1][iPix]));
occlVal+=normVarWeights[chn]*(thisDiff/pairVar/zMapVarThresholds[chn])/(pairDiff+k); // consider both - with /pairVar and without
pairVarDbg+=normVarWeights[chn]*pairVar;
pairDiffDbg+=normVarWeights[chn]*pairDiff;
thisDiffDbg+=normVarWeights[chn]*thisDiff;
double dd=slices[sIndex0][chn+1][1][iPix]-slices[sIndex1][chn+1][1][iPix];
thisDiffVar+=normVarWeights[chn]*Math.sqrt(0.5*(dd*dd+pairVar*pairVar));
// other variant - if pairDiff>1.0 - => occlVal=0.0;
}
if ((nPair==0) || (occlVarVal4) {
System.out.println( " pairVarDbgBest="+IJ.d2s(pairVarDbgBest,5));
}
}
aux0[i]=(float) occlVarVal;
// debug
// aux2[i]=(float) disparity;
aux3[i]=(float) pairVarDbgBest;
aux4[i]=(float) pairDiffDbgBest;
aux5[i]=(float) thisDiffDbgBest;
aux6[i]=(float) thisDiffVarBest;
}
zTile.setAux(0,aux0);
// debug
// zTile.setAux(2,aux2);
zTile.setAux(3,aux3);
zTile.setAux(4,aux4);
zTile.setAux(5,aux5);
zTile.setAux(6,aux6);
if (debugLevel>3){
float [][] debugData={aux2,aux3,aux4,aux5,aux6,aux0};
String [] debugTitles={"var","pair","this","occl","dvar","disp"};
(new ShowDoubleFloatArrays()).showArrays(
debugData,
paddedSize,
paddedSize,
true,
"R"+nImg+"_tile"+tileX+"-"+tileY,
debugTitles);
}
}
if (zMapCorrMask!=0) { // 8x8 correlation to find occlusion
double [][][]refineSmoothThis= new double [sImgSet.length][][];
double [][][]refineSmoothOther=new double [otherPairs.length][][];
for (int sIndex=0;sIndex3)?("RST"+nImg+"-"+sImg+"_"+chn):null);
}
}
for (int nPair=0;nPair3)?("RS-"+nImg+"-"+sImg0+"-"+sImg1+"_"+chn):null);
}
}
float [] aux1=new float[paddedLength];
//
float [] aux7=new float[paddedLength];
float [] aux8=new float[paddedLength];
float [] aux9=new float[paddedLength];
float [] aux10=new float[paddedLength];
for (int i=0;ibestOcclCorrVal)) {
bestOcclCorrVal=thisOcclCorrVal;
bestOcclCorrOther=refineSmoothOther[nPair][chn][iPix];
bestOcclCorrFirst=refineSmoothThis[sIndex0][chn][iPix];
bestOcclCorrSecond=refineSmoothThis[sIndex1][chn][iPix];
// bestPair=nPair;
}
}
occlCorrVal+=normCorrWeights[chn]*bestOcclCorrVal;
occlCorrOther+=normCorrWeights[chn]*bestOcclCorrOther;
occlCorrFirst+=normCorrWeights[chn]*bestOcclCorrFirst;
occlCorrSecond+=normCorrWeights[chn]*bestOcclCorrSecond;
}
// Save result to aux[1]
aux1[i]=(float) occlCorrVal;
aux7[i]=(float) occlCorrOther;
aux8[i]=(float) (occlCorrFirst+occlCorrSecond)/2;
aux9[i]=(float) occlCorrFirst;
aux10[i]=(float) occlCorrSecond;
}
zTile.setAux(1,aux1);
zTile.setAux(7,aux7);
zTile.setAux(8,aux8);
zTile.setAux(9,aux9);
zTile.setAux(10,aux10);
if (debugLevel>3){
float [][] debugData={aux1,aux7,aux8,aux9,aux10};
String [] debugTitles={"rslt","other","average","first","second"};
(new ShowDoubleFloatArrays()).showArrays(
debugData,
paddedSize,
paddedSize,
true,
"C"+nImg+"_tile"+tileX+"-"+tileY,
debugTitles);
}
}
return; // zTile.isForegroundValid(); // this current plane is not empty
}
public double [] normalizeWeights(int mask,double [] weights){
double [] maskedWeights=weights.clone();
for (int i=0;i2) System.out.println("setupZMap() zMapWOI.x="+zMapWOI.x+" zMapWOI.y="+zMapWOI.y+
" zMapWOI.width="+zMapWOI.width+" zMapWOI.height="+zMapWOI.height);
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger tileIndexAtomic = new AtomicInteger(0);
if (showProgress) IJ.showProgress(0.0);
if (showProgress) IJ.showStatus("Filtering zMap foreground for image "+nImg+" ("+(nImg+1)+" of "+this.disparityScales.length+") ...");
final int debugTile=debugRow*tilesX+debugColumn;
//TODO: define here
final double [] window=new double [4*this.overlapStep*this.overlapStep];
window[0]=Double.NaN;
final double [] refineWindow=new double [refineTilePeriod*refineTilePeriod*4];
refineWindow[0]=Double.NaN;
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
@Override
public void run() {
DoubleFHT doubleFHT = new DoubleFHT();
DoubleFHT refineFHT= new DoubleFHT();
for (int tile=tileIndexAtomic.getAndIncrement(); tile=debugThreshold) && (tile==debugTile);
int thisDebugLevel=debugLevel+(debugThis?2:0);
while (!filterForegroundZMapTile (
zMapWOI.x+tile%tilesX, //int tileX,
zMapWOI.y+tile/tilesX, //int tileY,
thisZMap,
nImg,
sImg, // list of second images
imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
imageFullWidth,
blurVarianceSigma,
refineTilePeriod,
refinePhaseCoeff,
refineHighPassSigma,
refineLowPassSigma,
refineCorrMaxDistance,
refineCorrThreshold,
refineSubPixel,
zMapMinForeground,
zMapVarMask,
zMapVarThresholds,
zMapVarWeights,
auxVarMode,
normalizeCorrelation,
zMapCorrMask,
zMapCorrThresholds,
zMapCorrWeights,
window,
doubleFHT,
refineWindow,
refineFHT,
threadsMax,
showProgress,
thisDebugLevel));
if (showProgress){
final int finalTile=tile;
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
IJ.showProgress(finalTile,tiles);
}
});
}
}
}
};
}
startAndJoin(threads);
IJ.showProgress(1.0);
}
private boolean filterForegroundZMapTile ( // false if nothing left in the current foreground, may repeat
int tileX,
int tileY,
ZTile [][] thisZMap,
int nImg,
int [] sImgSet, // list of second images
double [][][] imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
int imageFullWidth,
double blurVarianceSigma,
int refineTilePeriod,
double refinePhaseCoeff,
double refineHighPassSigma,
double refineLowPassSigma,
double refineCorrMaxDistance,
double refineCorrThreshold,
int refineSubPixel,
double zMapMinForeground,
int zMapVarMask,
double [] zMapVarThresholds,
double [] zMapVarWeights,
int auxVarMode,
int normalize,
int zMapCorrMask,
double [] zMapCorrThresholds,
double [] zMapCorrWeights,
double [] window,
DoubleFHT doubleFHT,
double [] refineWindow,
DoubleFHT refineFHT,
int threadsMax,
boolean showProgress,
int debugLevel){
// double [] normVarWeights= normalizeWeights(zMapVarMask,zMapVarWeights);
// double [] normCorrWeights=normalizeWeights(zMapCorrMask,zMapCorrWeights);
int auxChannelNumber=3;
/// final double [][] disparityScales=this.disparityScales;
// this filter will iterate through second images and AND filter for foreground pixels only
ZTile zTile=thisZMap[tileY][tileX];
if (zTile==null) return true; // tile is not defined, nothing left to do
// is there any forground?
if (!zTile.advanceForeground()) return true; // got to the end
if (zTile.enabledPixels[zTile.foregroundIndex]==null) {
zTile.initEnabledForegroundPixelsPlane();
}
int tileOverlap=zTile.getOverlap();
int paddedSize=zTile.getPaddedSize();
int paddedLength=zTile.getPaddedLength();
double disparity=zTile.maximums[zTile.foregroundIndex][0];
int size=this.overlapStep*2;
int margin=size/4-tileOverlap;
int [] dirs1={1,size+1,size,size-1,-1,-size-1,-size,-size+1,0};
for (int sIndex=0;sIndex3){ // +2 for selected tile
System.out.println("Debugging tileX="+tileX+", tileY="+tileY+" disparity="+disparity+" dX="+ dX+" dY="+dY);
int numActive=0;
for (int i=0;i=0;
df=Math.abs(df);
if ((d==0) || (df0.0){
(new DoubleGaussianBlur()) .blurDouble(
variance,
size,
size,
blurVarianceSigma,
blurVarianceSigma,
0.01);
}
for (int i=0;i=bDiff.length) || (i1>=variance.length)){
System.out.println("\nfilterForegroundZMapTile() i="+i+" margin="+margin+" paddedSize="+paddedSize+" bDiff.length="+bDiff.length+
" i1="+i1+" size="+size+" variance.length="+variance.length);
}
if ((bDiff[i1]/variance[i1])>zMapVarThresholds[chn]) { // 1024
zTile.disableForegroundPixel(i);
}
}
if (debugLevel>3){ // +2 for selected tile
String [] channelNames={"alpha","Y","Cb","Cr","Aux"};
String [] debugTitles={"bDiff","var","bdif/var","good","cumul"};
double [][] debugData=new double[5][length];
debugData[0]=bDiff;
debugData[1]=variance;
for (int i=0;izMapVarThresholds[chn])?-1.0:1.0;
debugData[4][i]=zTile.isEnabledForegroundPixel(i)?-1.0:1.0;
}
(new ShowDoubleFloatArrays()).showArrays(
debugData,
size,
size,
true,
"diff_"+channelNames[chn+1]+"_"+tileX+"-"+tileY+"_"+nImg+"-"+sImg,
debugTitles);
}
} else { // aux channel
for (int i=0;i3){ // +2 for selected tile
(new ShowDoubleFloatArrays()).showArrays(
bDiff,
size,
size,
"AUX_"+nImg+"-"+sImg+"__"+tileX+"-"+tileY);
}
for (int i=0;i3)?("CORR"+nImg+"-"+sImg+"_"+chn):null);
for (int i=0;i2) System.out.println("setupZMap() zMapWOI.x="+zMapWOI.x+" zMapWOI.y="+zMapWOI.y+
" zMapWOI.width="+zMapWOI.width+" zMapWOI.height="+zMapWOI.height);
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger tileIndexAtomic = new AtomicInteger(0);
if (showProgress) IJ.showProgress(0.0);
if (showProgress) IJ.showStatus("refinePlaneDisparity for image "+nImg+" ("+(nImg+1)+" of "+this.disparityScales.length+") ...");
final int debugTile=debugRow*tilesX+debugColumn;
//TODO: define here
// final double [] window=new double [4*this.overlapStep*this.overlapStep];
// window[0]=Double.NaN;
// final double [] refineWindow=new double [refineTilePeriod*refineTilePeriod*4];
// refineWindow[0]=Double.NaN;
final int size=2*this.overlapStep;
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
@Override
public void run() {
DoubleFHT doubleFHT = new DoubleFHT();
final double [] window=new double [size*size];
window[0]=Double.NaN;
// DoubleFHT refineFHT= new DoubleFHT();
for (int tile=tileIndexAtomic.getAndIncrement(); tile=debugThreshold) && (tile==debugTile);
int thisDebugLevel=debugLevel+(debugThis?2:0);
refinePlaneDisparityTile (
zMapWOI.x+tile%tilesX, //int tileX,
zMapWOI.y+tile/tilesX, //int tileY,
zMapFinal,
nImg,
sImg, // list of second images
imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
imageFullWidth,
// blurVarianceSigma,
// refineTilePeriod,
refinePhaseCoeff,
refineHighPassSigma,
refineLowPassSigma,
refineCorrMaxDistance,
refineCorrThreshold,
refineSubPixel,
// zMapMinForeground,
// zMapVarMask,
// zMapVarThresholds,
// zMapVarWeights,
// auxVarMode,
// normalizeCorrelation,
zMapCorrMask,
zMapCorrThresholds,
zMapCorrWeights,
window,
doubleFHT,
// refineWindow,
// refineFHT,
disparityMax,
disparityMin,
minAbsolute, // or NaN - will use enabled/disabled state of the tile
minRelative,
filterByForeground, // apply known certain masks
filterByForegroundMargin,
filterByDisabled,
disparityTolearnce,
maskBlurSigma,
corrHighPassSigma, // subtract blurred version to minimize correlation caused by masks
threadsMax,
showProgress,
thisDebugLevel);
if (showProgress){
final int finalTile=tile;
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
IJ.showProgress(finalTile,tiles);
}
});
}
}
}
};
}
startAndJoin(threads);
IJ.showProgress(1.0);
}
private void refinePlaneDisparityTile ( // false if nothing left in the current foreground, may repeat
int tileX,
int tileY,
ZTile [][][] allZMap,
int nImg,
int [] sImgSet, // list of second images
double [][][] imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
int imageFullWidth,
// double blurVarianceSigma,
// int refineTilePeriod,
double refinePhaseCoeff,
double refineHighPassSigma,
double refineLowPassSigma,
double refineCorrMaxDistance,
double refineCorrThreshold,
int refineSubPixel,
// double zMapMinForeground,
// int zMapVarMask,
// double [] zMapVarThresholds,
// double [] zMapVarWeights,
// int auxVarMode,
// int normalize,
int zMapCorrMask,
double [] zMapCorrThresholds,
double [] zMapCorrWeights,
double [] window,
DoubleFHT doubleFHT,
// double [] refineWindow,
// DoubleFHT refineFHT,
// new arguments
// int combineMode, // different image pairs - 0
double disparityMax,
double disparityMin,
double minAbsolute, // or NaN - will use enabled/disabled state of the tile
double minRelative,
boolean filterByForeground, // apply known certain masks
double filterByForegroundMargin,
boolean filterByDisabled,
double disparityTolearnce,
double maskBlurSigma,
double corrHighPassSigma, // subtract blurred version to minimize correlation caused by masks
int threadsMax,
boolean showProgress,
int debugLevel){
// TODO: also calculate "unlikely" - high autocorrelation, not occluded, low inter-correlation
double [] normCorrWeights=normalizeWeights(zMapCorrMask,zMapCorrWeights);
ZTile zTile=allZMap[nImg][tileY][tileX];
int tileOverlap=zTile.getOverlap();
int paddedSize=zTile.getPaddedSize();
int paddedLength=zTile.getPaddedLength();
int size=this.overlapStep*2;
int length=size*size;
double[] zeros=new double [length];
for (int i=0;i3) {
System.out.println ("zTile.setMinCorrelations("+minAbsolute+","+ minRelative+")");
boolean [] enabled=zTile.enabledPlane;
for (int plane=0;plane3) {
System.out.println ("refinePlaneDisparityTile("+tileX+","+tileY+", ...) plane="+plane+
" disparity="+zTile.getPlaneDisparity(plane)+" strength="+zTile.getPlaneStrength(plane));
}
if (!(zTile.getPlaneDisparity(plane)disparityMax)){ // NaN is OK
double disparity= zTile.getPlaneDisparity(plane);
if (debugLevel>3) {
// double testDisparityError=0.5;
System.out.println (" processing plane "+plane+" nImg="+nImg+
" disparity="+disparity+
" filterByForegroundMargin="+filterByForegroundMargin);
// disparity+=testDisparityError;
// System.out.println ("**** modified disparity for testing, new disparity="+disparity);
}
double [] thisEnabledMask;
if (filterByForeground){
boolean [] bMask=zTile.getEnabledNonOccluded(
disparity+filterByForegroundMargin,
disparity,
filterByDisabled?disparityTolearnce:Double.NaN,
debugLevel);
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
bMask,
paddedSize,
paddedSize,
"bMask_N"+nImg+"-X"+tileX+"-Y"+tileY
);
}
thisEnabledMask=zeros.clone();
for (int i=0;i0.0){
(new DoubleGaussianBlur()).blurDouble(
thisEnabledMask,
size,
size,
maskBlurSigma,
maskBlurSigma,
0.01);
}
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
thisEnabledMask,
size,
size,
"thisEnabledMask_N"+nImg+"-X"+tileX+"-Y"+tileY
);
}
} else {
thisEnabledMask=null;
}
double [][] pairEnabledMask=new double [sImgSet.length][];
double [] disparityArray=null;
// int disparityRange=(int) Math.ceil(refineCorrMaxDistance*refineSubPixel*2);
for (int sIndex=0;sIndex3){ // +2 for selected tile
System.out.println("Iterating tileX="+tileX+", tileY="+tileY+" disparity="+disparity+" dX="+ dX+" dY="+dY+
" sIndex="+sIndex+" nImg="+nImg+" sImg="+sImg);
}
if (filterByForeground){
double [] otherNonOccluded= getNonOccluded(
allZMap[sImg],//ZTile [][] thisZMap,
tileX*this.overlapStep+this.overlapStep/2+dX, //double xc,
tileY*this.overlapStep+this.overlapStep/2+dY, //double yc,
disparity, //double disparity,
filterByForegroundMargin,
filterByDisabled?disparityTolearnce:Double.NaN,
debugLevel);
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
otherNonOccluded,
paddedSize,
paddedSize,
"PEM_N"+nImg+"-S"+sImg+"-X"+tileX+"-Y"+tileY
);
}
pairEnabledMask[sIndex]=zeros.clone();
for (int i=0;i0.0){
(new DoubleGaussianBlur()).blurDouble(
pairEnabledMask[sIndex],
size,
size,
maskBlurSigma,
maskBlurSigma,
0.01);
}
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
pairEnabledMask[sIndex],
size,
size,
"pairEnabledMask_N"+nImg+"-S"+sImg+"-X"+tileX+"-Y"+tileY
);
}
//multiply masks from both images
for (int i=0;i3){ // +2 for selected tile
System.out.println("Debugging tileX="+tileX+", tileY="+tileY+" plane="+plane+" disparity="+disparity+" dX="+ dX+" dY="+dY+
" sIndex="+sIndex+" nImg="+nImg+" sImg="+sImg);
int numActive=0;
for (int i=0;i0.0){
double [] loPass=pairs[chn][n].clone();
(new DoubleGaussianBlur()).blurDouble(
loPass,
size,
size,
corrHighPassSigma,
corrHighPassSigma,
0.01);
for (int i=0;i3){ // +2 for selected tile
System.out.println("centerCorr[0]="+centerCorr[0]+" centerCorr[1]="+centerCorr[1]+" centerCorr[2]="+centerCorr[2]);
// show high-pass/masked pairs
int numActive=0;
for (int i=0;i=debugData.length) || ((i-1)>=pairs.length)){
System.out.println(
" index="+index+
" debugData.length="+debugData.length+
" i="+i+
" pairs.length="+pairs.length+
" slices.length="+slices.length+
" slices[sIndex].length="+slices[sIndex].length
);
}
debugData[index++]=pairs[i-1][0]; //ava.lang.ArrayIndexOutOfBoundsException: 3
debugTitles[index]=channelNames[i]+"-"+sImg;
debugData[index++]=pairs[i-1][1];
}
(new ShowDoubleFloatArrays()).showArrays(
debugData,
size,
size,
true,
"P"+plane+"_N"+nImg+"-"+sImg+"_sIndex-"+sIndex+"_tile"+tileX+"-"+tileY,
debugTitles);
}
// TODO: now look if the correlation is strong enough to perform correction (minimal of 3?)
double minCorr=1.0;
for (int i=0;i3){ // +2 for selected tile
System.out.println("centerCorr[0]="+centerCorr[0]+" centerCorr[1]="+centerCorr[1]+" centerCorr[2]="+centerCorr[2]+
" minCorr="+minCorr+" refineCorrThreshold="+refineCorrThreshold);
}
if (minCorr>=refineCorrThreshold){
double [] combinedChnCorr=new double [length];
for (int i=0;imax){
max=combinedChnCorr[i];
iMax=i;
}
int ixc=iMax%size-size/2;
int iyc=iMax/size-size/2;
if (debugLevel>3) System.out.println("refineCorrelation(): max="+max+" iMax="+iMax+" refineCorrMaxDistance="+refineCorrMaxDistance+" corrMaxDist2="+corrMaxDist2+
" r2="+(ixc*ixc+iyc*iyc)+" r="+Math.sqrt(ixc*ixc+iyc*iyc));
if ((ixc*ixc+iyc*iyc)<=corrMaxDist2){ // maximum close enough
double [] upsampled=doubleFHT.upsample(combinedChnCorr,refineSubPixel);
int interpolatedSize=size*refineSubPixel;
int interpolatedCenter=(interpolatedSize+1)*interpolatedSize/2;
if (debugLevel>3) {
(new ShowDoubleFloatArrays()).showArrays(
upsampled,
interpolatedSize,
interpolatedSize,
"U"+plane+"_N"+nImg+"-"+sImg+"_sIndex-"+sIndex+"_tile"+tileX+"-"+tileY);
}
if (disparityArray==null) {
disparityArray=new double [2*halfDispRange+1];
for (int i=0;i5){
System.out.println("disparityArray["+i+"]="+disparityArray[i]+
" deltaX="+deltaX+" deltaY="+deltaY+" iX="+iX+" iY="+iY+" index00="+index00+" index01="+index01+
" index10="+index10+" index11="+index11);
}
}
}
} // if (minCorr>=refineCorrThreshold)
} //for (int sIndex=0;sIndexdisparityArray[iMax]) iMax=i;
}
double disparCorr=((double) (iMax-halfDispRange))/refineSubPixel;
if (debugLevel>5){
for (int i=0;i3){
System.out.println("Old disparity["+plane+"]="+disparity+" new disparity="+(disparity+disparCorr));
}
zTile.setPlaneDisparity(disparity+disparCorr,plane);
}
}
}
} // end for (int plane...)
}
public void planeLikely(
// int nImg,
// int [] sImg, // list of second images
double [][][] imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
int imageFullWidth,
double blurVarianceSigma,
int refineTilePeriod,
double subTilePhaseCoeff,
double subTileHighPassSigma,
double subTileLowPassSigma,
// double refineCorrMaxDistance,
// double refineCorrThreshold,
int refineSubPixel,
double zMapMinForeground,
int zMapVarMask,
double [] zMapVarThresholds,
double [] zMapVarWeights,
int auxVarMode,
int normalizeCorrelation,
int zMapCorrMask,
double [] zMapCorrThresholds,
double [] zMapCorrThresholdsRel,
double [] zMapCorrWeights,
int debugRow,
int debugColumn,
int combineMode, // different image pairs - 0
double disparityMax,
double disparityMin,
double minAbsolute, // or NaN - will use enabled/disabled state of the tile
double minRelative,
boolean filterByForeground, // apply known certain masks
double filterByForegroundMargin,
boolean filterByDisabled,
double disparityTolearnce,
double maskBlurSigma,
double corrHighPassSigma, // subtract blurred version to minimize correlation caused by masks
double varianceBlurScale,
double kLocal,
int matchStatMode,
int threadsMax,
boolean showProgress,
int debugLevel){
// process all image pairs
for (int nImg=0;nImg2) System.out.println("setupZMap() zMapWOI.x="+zMapWOI.x+" zMapWOI.y="+zMapWOI.y+
" zMapWOI.width="+zMapWOI.width+" zMapWOI.height="+zMapWOI.height);
final Thread[] threads = newThreadArray(threadsMax);
final AtomicInteger tileIndexAtomic = new AtomicInteger(0);
if (showProgress) IJ.showProgress(0.0);
if (showProgress) IJ.showStatus("planeLikely for image "+nImg+" ("+(nImg+1)+" of "+this.disparityScales.length+") ...");
final int debugTile=debugRow*tilesX+debugColumn;
//TODO: define here
// final double [] window=new double [4*this.overlapStep*this.overlapStep];
// window[0]=Double.NaN;
// final double [] refineWindow=new double [refineTilePeriod*refineTilePeriod*4];
// refineWindow[0]=Double.NaN;
final int size=2*this.overlapStep;
for (int ithread = 0; ithread < threads.length; ithread++) {
threads[ithread] = new Thread() {
@Override
public void run() {
DoubleFHT doubleFHT = new DoubleFHT();
final double [] window=new double [size*size];
window[0]=Double.NaN;
final double [] subTileWindow=new double [refineTilePeriod*refineTilePeriod*4];
subTileWindow[0]=Double.NaN;
DoubleFHT subTileFHT= new DoubleFHT();
for (int tile=tileIndexAtomic.getAndIncrement(); tile=debugThreshold) && (tile==debugTile);
int thisDebugLevel=debugLevel+(debugThis?2:0);
planeLikelyTile (
zMapWOI.x+tile%tilesX, //int tileX,
zMapWOI.y+tile/tilesX, //int tileY,
zMapFinal,
nImg,
sImg, // list of second images
imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
imageFullWidth,
blurVarianceSigma,
// refineTilePeriod,
subTilePhaseCoeff,
subTileHighPassSigma,
subTileLowPassSigma,
// refineCorrMaxDistance,
// refineCorrThreshold,
refineSubPixel,
zMapMinForeground,
zMapVarMask,
zMapVarThresholds,
zMapVarWeights,
auxVarMode,
normalizeCorrelation,
zMapCorrMask,
zMapCorrThresholds,
zMapCorrThresholdsRel,
zMapCorrWeights,
window,
doubleFHT,
subTileWindow,
subTileFHT,
combineMode, // different image pairs - 0
disparityMax,
disparityMin,
minAbsolute, // or NaN - will use enabled/disabled state of the tile
minRelative,
filterByForeground, // apply known certain masks
filterByForegroundMargin,
filterByDisabled,
disparityTolearnce,
maskBlurSigma,
corrHighPassSigma, // subtract blurred version to minimize correlation caused by masks
varianceBlurScale,
kLocal,
matchStatMode,
threadsMax,
showProgress,
thisDebugLevel);
if (showProgress){
final int finalTile=tile;
SwingUtilities.invokeLater(new Runnable() {
@Override
public void run() {
IJ.showProgress(finalTile,tiles);
}
});
}
}
}
};
}
startAndJoin(threads);
IJ.showProgress(1.0);
}
private void planeLikelyTile ( // false if nothing left in the current foreground, may repeat
int tileX,
int tileY,
ZTile [][][] allZMap,
int nImg,
int [] sImgSet, // list of second images
double [][][] imageData, // [img]{Alpha, Y,Cb,Cr, Ext}. Alpha may be null - will handle later?
int imageFullWidth,
double blurVarianceSigma,
// int refineTilePeriod,
double subTilePhaseCoeff,
double subTileHighPassSigma,
double subTileLowPassSigma,
// double refineCorrMaxDistance,
// double refineCorrThreshold,
int refineSubPixel,
double zMapMinForeground,
int zMapVarMask,
double [] zMapVarThresholds,
double [] zMapVarWeights,
int auxVarMode,
int normalize,
int zMapCorrMask,
double [] zMapCorrThresholds,
double [] zMapCorrThresholdsRel,
double [] zMapCorrWeights,
double [] window,
DoubleFHT doubleFHT,
double [] subTileWindow,
DoubleFHT subTileFHT,
// new arguments
int combineMode, // different image pairs - 0
double disparityMax,
double disparityMin,
double minAbsolute, // or NaN - will use enabled/disabled state of the tile
double minRelative,
boolean filterByForeground, // apply known certain masks
double filterByForegroundMargin,
boolean filterByDisabled,
double disparityTolearnce,
double maskBlurSigma,
double corrHighPassSigma, // subtract blurred version to minimize correlation caused by masks
//tone-matching statistical parameters
double varianceBlurScale,
double kLocal,
int matchStatMode,
int threadsMax,
boolean showProgress,
int debugLevel){
// TODO: also calculate "unlikely" - high autocorrelation, not occluded, low inter-correlation
double [] normVarWeights= normalizeWeights(zMapVarMask,zMapVarWeights);
double [] normCorrWeights=normalizeWeights(zMapCorrMask,zMapCorrWeights);
// int auxChannelNumber=3;
ZTile zTile=allZMap[nImg][tileY][tileX];
int tileOverlap=zTile.getOverlap();
int paddedSize=zTile.getPaddedSize();
int paddedLength=zTile.getPaddedLength();
int size=this.overlapStep*2;
int length=size*size;
double[] zeros=new double [length];
for (int i=0;i3) {
System.out.println ("zTile.setMinCorrelations("+minAbsolute+","+ minRelative+")");
boolean [] enabled=zTile.enabledPlane;
for (int plane=0;plane3) {
System.out.println ("planeLikelyTile("+tileX+","+tileY+", ...) plane="+plane+
" disparity="+zTile.getPlaneDisparity(plane)+" strength="+zTile.getPlaneStrength(plane));
}
if (!(zTile.getPlaneDisparity(plane)disparityMax)){ // NaN is OK
planeStrength[plane]=new double [sImgSet.length][];
double disparity= zTile.getPlaneDisparity(plane);
if (debugLevel>3) {
System.out.println (" processing plane "+plane+" nImg="+nImg+
" disparity="+disparity+
" filterByForegroundMargin="+filterByForegroundMargin);
}
double [] thisEnabledMask;
if (filterByForeground){
boolean [] bMask=zTile.getEnabledNonOccluded(
disparity+filterByForegroundMargin,
disparity,
filterByDisabled?disparityTolearnce:Double.NaN,
debugLevel);
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
bMask,
paddedSize,
paddedSize,
"bMask_N"+nImg+"-X"+tileX+"-Y"+tileY
);
}
thisEnabledMask=zeros.clone();
for (int i=0;i0.0){
(new DoubleGaussianBlur()).blurDouble(
thisEnabledMask,
size,
size,
maskBlurSigma,
maskBlurSigma,
0.01);
}
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
thisEnabledMask,
size,
size,
"thisEnabledMask_N"+nImg+"-X"+tileX+"-Y"+tileY
);
}
} else {
thisEnabledMask=null;
}
double [][] pairEnabledMask=new double [sImgSet.length][];
// double [] disparityArray=null; // needed??
for (int sIndex=0;sIndex3){ // +2 for selected tile
System.out.println("Iterating tileX="+tileX+", tileY="+tileY+" disparity="+disparity+" dX="+ dX+" dY="+dY+
" sIndex="+sIndex+" nImg="+nImg+" sImg="+sImg);
}
if (filterByForeground){
double [] otherNonOccluded= getNonOccluded(
allZMap[sImg],//ZTile [][] thisZMap,
tileX*this.overlapStep+this.overlapStep/2+dX, //double xc,
tileY*this.overlapStep+this.overlapStep/2+dY, //double yc,
disparity, //double disparity,
filterByForegroundMargin,
filterByDisabled?disparityTolearnce:Double.NaN,
debugLevel);
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
otherNonOccluded,
paddedSize,
paddedSize,
"PEM_N"+nImg+"-S"+sImg+"-X"+tileX+"-Y"+tileY
);
}
pairEnabledMask[sIndex]=zeros.clone();
for (int i=0;i0.0){
(new DoubleGaussianBlur()).blurDouble(
pairEnabledMask[sIndex],
size,
size,
maskBlurSigma,
maskBlurSigma,
0.01);
}
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
pairEnabledMask[sIndex],
size,
size,
"pairEnabledMask_N"+nImg+"-S"+sImg+"-X"+tileX+"-Y"+tileY
);
}
//multiply masks from both images
for (int i=0;i3){ // +2 for selected tile
System.out.println("Debugging tileX="+tileX+", tileY="+tileY+" plane="+plane+" disparity="+disparity+" dX="+ dX+" dY="+dY+
" sIndex="+sIndex+" nImg="+nImg+" sImg="+sImg);
int numActive=0;
for (int i=0;i0.0){
double [] loPass=pairs[chn][n].clone();
(new DoubleGaussianBlur()).blurDouble(
loPass,
size,
size,
corrHighPassSigma,
corrHighPassSigma,
0.01);
for (int i=0;i3){ // +2 for selected tile
// System.out.println("centerCorr[0]="+centerCorr[0]+" centerCorr[1]="+centerCorr[1]+" centerCorr[2]="+centerCorr[2]);
// show high-pass/masked pairs
int numActive=0;
for (int i=0;i=debugData.length) || ((i-1)>=pairs.length)){
System.out.println(
" index="+index+
" debugData.length="+debugData.length+
" i="+i+
" pairs.length="+pairs.length+
" slices.length="+slices.length+
" slices[sIndex].length="+slices[sIndex].length
);
}
debugData[index++]=pairs[i-1][0]; //ava.lang.ArrayIndexOutOfBoundsException: 3 //Exception in thread "Thread-1033" java.lang.NullPointerException
debugTitles[index]=channelNames[i]+"-"+sImg;
debugData[index++]=pairs[i-1][1];
}
(new ShowDoubleFloatArrays()).showArrays(
debugData,
size,
size,
true,
"P"+plane+"_N"+nImg+"-"+sImg+"_sIndex-"+sIndex+"_tile"+tileX+"-"+tileY,
debugTitles);
}
//=================== several alternative modes. For likelyhood - find linear relation between both images using weights - proportional to cross-correlation
double [] linearMatchWeight=pairEnabledMask[sIndex].clone();
if (corrs!=null){ // and some mode?
int tileSize=(int) Math.sqrt(subTileWindow.length); // 8
int tileStep=tileSize/2; //4
int tileExtra=size-tileSize; // 24??
// int tileLength=tileSize*tileSize; // 64
// int tileCenter=(tileSize+1)*tileSize/2; // 36
int topLeft=(size+1)*(size/4-tileSize/4);
int numTilesRow=size/(2*tileStep); // 4
// accumulate result
double [] corrTile=new double [length];
for (int i=0;i3){
System.out.println ("corrs["+subTileY+"]["+subTileX+"]={"+
corrs[subTileY][subTileX][0]+","+
corrs[subTileY][subTileX][1]+","+
corrs[subTileY][subTileX][2]+"} aCorr="+aCorr+
" corr="+corr);
}
for (int iy=0;iy0.0)?corrTile[i]:0.0);
double sumWindow=0.0;
double sumLMW=0.0;
for (int i=0;i3){ // +2 for selected tile
double [][] debugData={linearMatchWeight,linearMatchWeightNorm};
String [] debugTitles={"LMW","Norm"};
(new ShowDoubleFloatArrays()).showArrays(
debugData,
size,
size,
true,
"LMW"+plane+"_N"+nImg+"-"+sImg+"_sIndex-"+sIndex+"_tile"+tileX+"-"+tileY,
debugTitles);
}
double [][]staging=new double [slices[sIndex].length-1][];
double [][][]variance=new double[slices[sIndex].length-1][][];
double [] avgVariance=new double[slices[sIndex].length-1];
double [] avgRange= new double[slices[sIndex].length-1];
for (int chn=0;chn5); //double var2){ // second image variance
}
double [] fVar=new double[2];
for (int n=0;n<2;n++){
fVar[n]=S2tot[n]/S0tot[n]-S1tot[n]*S1tot[n]/S0tot[n]/S0tot[n];
}
avgVariance[chn]=Math.sqrt(sumWV2/sumW);
avgRange[chn]= Math.sqrt(0.5*(fVar[0]+fVar[1]));
if (debugLevel>3){
System.out.println("Average variance for plane "+plane+" nImg="+nImg+" sImg="+sImg+
" chn="+chn+" is "+avgVariance[chn]+" total tile variance="+avgRange[chn]);
}
// blur stage with scaled variance
this.photometric.blurStaging(
staging[chn],
chn,
avgVariance[chn]*varianceBlurScale);
//varianceBlurScale
}else{
staging[chn]=null;
variance[chn]=null;
avgVariance[chn]=0.0;
avgRange[chn]=0.0;
}
}
double sumRelVar=0.0;
for (int chn=0;chn0.0){
if (debugLevel>3){
System.out.println("avgVariance["+chn+"]="+avgVariance[chn]+
" avgRange["+chn+"]="+avgRange[chn]+
" getAverageVariance("+nImg+","+chn+")="+this.photometric.getAverageVariance(nImg,chn));
}
// avgVariance[chn]/=this.photometric.getAverageVariance(nImg,chn);
// avgVariance[chn]*=avgVariance[chn];
avgRange[chn]/=this.photometric.getAverageVariance(nImg,chn);
avgRange[chn]*=avgRange[chn];
sumRelVar+=avgVariance[chn];
}
for (int chn=0;chn0.0){
avgVariance[chn]/=sumRelVar;
if (debugLevel>3){
System.out.println(" normalized avgVariance["+chn+"]="+avgVariance[chn]);
}
}
if (debugLevel>3){
int numActive=0;
for (int i=0;i3) planeStrengthChn=new double [staging.length][length];
for (int chn=0;chn3) {
planeStrengthChn[chn]=new double[length];
for (int i=0;i0){
double strength=this.photometric.getStrength(
staging[chn],
chn,
slices[sIndex][chn+1][0][i], //double v1,
slices[sIndex][chn+1][1][i]); //double v2
strength*=linearMatchWeight[i]; // remove??
if (window[i]>0.0) strength/=window[i];
planeStrength[plane][sIndex][i]+=avgVariance[chn]*normVarWeights[chn]*strength;
if (debugLevel>3) planeStrengthChn[chn][i]=strength;
}
} else {
if (debugLevel>3) planeStrengthChn[chn]=null;
}
}
if (debugLevel>3){
int numActive=0;
for (int i=0;i0.0) planeStrengthCombo[plane][i]+=planeStrength[plane][sIndex][i];
}
break;
case 2: // min
for (int i=0;i planeStrengthCombo[plane][i]) planeStrengthCombo[plane][i]=planeStrength[plane][sIndex][i];
}
break;
}
}
double a=1.0/planeStrength[plane].length;
switch (matchStatMode) {
case 0: // multiply
for (int i=0;i0.0) planeStrengthCombo[plane][i]=Math.pow(planeStrengthCombo[plane][i],a);
}
break;
case 1: // add
for (int i=0;i0) && // should be positive
((planeIndex[i]<0) ||((planeStrengthCombo[planeIndex[i]]!=null) &&
(planeStrengthCombo[plane]!=null) &&
(planeStrengthCombo[planeIndex[i]][iPix] < planeStrengthCombo[plane][iPix])))) { //java.lang.NullPointerException
planeIndex[i]=plane;
numUsed_dbg++;
}
}
if (debugLevel>3) {
System.out.println ("planeLikelyTile("+tileX+","+tileY+", ...) plane="+plane+
" disparity="+zTile.getPlaneDisparity(plane)+" strength="+zTile.getPlaneStrength(plane)+
" used "+numUsed_dbg+"pixels"); // only scan inner pixels?
}
}
if (debugLevel>3) {
int numActivePlanes=0;
for (int plane=0;plane0){
int index=0;
int [] planeNumbers=new int [numActivePlanes];
for (int plane=0;plane=0)?zTile.getPlaneDisparity( planeIndex[i]):0.0); // make it NaN?
zTile.setAux(0,aux0);
}
/**
* calculate auto/inter correaltion of the smaller subtiles
* @param doubleFHT - DHT instance reused by the thread (should be initialized)
* @param data first/second array - now [2][32*32]
* @param window - subtile window (defines subtile size) - now 8x8=64 - has to be initialized to an array with [0] equal to NaN
* @param phaseCoeff phase correlation fraction (0 - normal, 1.0 - pure phase)
* @param highPassSigma high-pass sigma - fraction of the full range
* @param lowPassSigma low-pass sigma, in pixels
* @param debugLevel
* @return for each row and column - value of the correlation in the center (auto 00, auto 11, inter 01): [vert][hor]{corr00, corr11, corr01}
*/
public double [][][] subTileCorrelation(
DoubleFHT doubleFHT, // will generate if null, can be used to reuse initialization
double [][] data, // difines size
double [] window, // defines tile size
double phaseCoeff, // phase correlation fraction (0 - normal, 1.0 - pure phase)
double highPassSigma,
double lowPassSigma,
int debugLevel){
int size=(int) Math.sqrt(data[0].length); // 32
int tileSize=(int) Math.sqrt(window.length); // 8
int tileStep=tileSize/2; //4
int tileExtra=size-tileSize; // 24??
int tileLength=tileSize*tileSize; // 64
// int tileCenter=(tileSize+1)*tileSize/2; // 36
int numTilesRow=size/(2*tileStep); // 4
double [][][] corrs=new double[numTilesRow][numTilesRow][3];
if (Double.isNaN(window[0])){
double [] window1d=doubleFHT.getHamming1d(tileSize); // Hamming
for (int i=0;iiXYtl[0]/this.overlapStep,
iXYbr[1]/this.overlapStep>iXYtl[1]/this.overlapStep};
int tileX=iXYtl[0]/this.overlapStep;
int tileY=iXYtl[1]/this.overlapStep;
boolean [] bothXY={
((int) Math.floor(((double) iXYbr[0])/this.overlapStep))>((int) Math.floor(((double)iXYtl[0])/this.overlapStep)),
((int) Math.floor(((double) iXYbr[1])/this.overlapStep))>((int) Math.floor(((double)iXYtl[1])/this.overlapStep))};
*/
boolean [] bothXY={
Math.floor(((double) iXYbr[0])/this.overlapStep)>Math.floor(((double)iXYtl[0])/this.overlapStep),
Math.floor(((double) iXYbr[1])/this.overlapStep)>Math.floor(((double)iXYtl[1])/this.overlapStep)};
int tileX=(int)Math.floor(((double) iXYtl[0])/this.overlapStep);
int tileY=(int)Math.floor(((double) iXYtl[1])/this.overlapStep);
if (debugLevel>3) { //(tileX<0) || (tileY<0)){
System.out.println("getNonOccluded(thisZMap,"+xc+","+yc+","+disparity+
","+foregroundDisparityMargin+","+disparityTolerance+","+debugLevel+
") :iXYtl={"+iXYtl[0]+","+iXYtl[1]+"} :iXYbr={"+iXYbr[0]+","+iXYbr[1]+"} :dXY={"+dXY[0]+","+dXY[1]+"} tileX="+
tileX+" tileY="+tileY);
}
double disparityThresholdFg=disparity+foregroundDisparityMargin;
// enMasks[0][0]=thisZMap[tileY][tileX].getNonOccluded(disparityThresholdFg);
// enMasks[0][1]=(bothXY[0])?thisZMap[tileY][tileX+1].getNonOccluded(disparityThresholdFg):null;
// enMasks[1][0]=(bothXY[1])?thisZMap[tileY+1][tileX].getNonOccluded(disparityThresholdFg):null;
// enMasks[1][0]=(bothXY[0]&&bothXY[1])?thisZMap[tileY+1][tileX+1].getNonOccluded(disparityThresholdFg):null;
enMasks[0][0]=getEnabledNonOccluded(
thisZMap,
tileX,
tileY,
disparityThresholdFg,
disparity,
disparityTolerance,
debugLevel);
enMasks[0][1]=(bothXY[0])?
getEnabledNonOccluded(
thisZMap,
tileX+1,
tileY,
disparityThresholdFg,
disparity,
disparityTolerance,
debugLevel):null;
enMasks[1][0]=(bothXY[1])?
getEnabledNonOccluded(
thisZMap,
tileX,
tileY+1,
disparityThresholdFg,
disparity,
disparityTolerance,
debugLevel):null;
enMasks[1][1]=(bothXY[0]&&bothXY[1])?getEnabledNonOccluded(
thisZMap,
tileX+1,
tileY+1,
disparityThresholdFg,
disparity,
disparityTolerance,
debugLevel):null;
// int paddedSize=thisZMap[tileY][tileX].getPaddedSize();
// int margin=this.overlapStep/2-thisZMap[tileY][tileX].getOverlap();
// int padding=thisZMap[tileY][tileX].getOverlap();
int padding=getPadding();
int [][] whereX=new int [this.paddedSize+(bothXY[0]?1:0)][2];
int [][] whereY=new int [this.paddedSize+(bothXY[1]?1:0)][2];
int [] tile00tl= {tileX*this.overlapStep-padding,tileY*this.overlapStep-padding};
for (int i=0;i=(tile00tl[0]+this.overlapStep))){
whereX[i][0]-=this.overlapStep;
whereX[i][1]=1;
}
}
for (int i=0;i=(tile00tl[1]+this.overlapStep))){
whereY[i][0]-=this.overlapStep;
whereY[i][1]=1;
}
}
int longPaddedSize=this.paddedSize+1;
boolean [] extraBits=new boolean[longPaddedSize*longPaddedSize];
for (int y=0;y=enMasks[whereY[y][1]][whereX[x][1]].length){
System.out.println(
" xc="+xc+
" yc="+yc+
" tileX="+tileX+
" tileY="+tileY+
" iXYtl={"+iXYtl[0]+","+iXYtl[1]+"}"+
" iXYbr={"+iXYbr[0]+","+iXYbr[1]+"}"+
" bothXY={"+bothXY[0]+","+bothXY[1]+"}"+
" whereX.length="+whereX.length+
" whereY.length="+whereY.length+
" whereY["+y+"][0]="+whereY[y][0]+
" whereY["+y+"][1]="+whereY[y][1]+
" whereX["+x+"][0]="+whereX[x][0]+
" whereX["+x+"][1]="+whereX[x][1]+
" enMasks["+whereY[y][1]+"]["+whereX[x][1]+"].length="+enMasks[whereY[y][1]][whereX[x][1]].length+
" whereY[y][0]*this.paddedSize + whereX[x][0]="+(whereY[y][0]*this.paddedSize + whereX[x][0])+
" Math.floor(((double)iXYtl[0])/this.overlapStep)="+Math.floor(((double)iXYtl[0])/this.overlapStep)+
" Math.floor(((double)iXYbr[0])/this.overlapStep)="+Math.floor(((double)iXYbr[0])/this.overlapStep)+
" Math.floor(((double)iXYtl[1])/this.overlapStep)="+Math.floor(((double)iXYtl[1])/this.overlapStep)+
" Math.floor(((double)iXYbr[1])/this.overlapStep)="+Math.floor(((double)iXYbr[1])/this.overlapStep)+
" (int) Math.floor(((double)iXYtl[0])/this.overlapStep)="+((int) Math.floor(((double)iXYtl[0])/this.overlapStep))+
" (int) Math.floor(((double)iXYbr[0])/this.overlapStep)="+((int) Math.floor(((double)iXYbr[0])/this.overlapStep))+
" (int) Math.floor(((double)iXYtl[1])/this.overlapStep)="+((int) Math.floor(((double)iXYtl[1])/this.overlapStep))+
" (int) Math.floor(((double)iXYbr[1])/this.overlapStep)="+((int) Math.floor(((double)iXYbr[1])/this.overlapStep))/*+
" ((int) (Math.floor((double)iXYbr[1])/this.overlapStep))>((int) (Math.floor((double)iXYtl[1])/this.overlapStep))="+
(((int) (Math.floor((double)iXYbr[1])/this.overlapStep))>((int) (Math.floor((double)iXYtl[1])/this.overlapStep)))+
" ((Math.floor((double)iXYbr[1])/this.overlapStep))>((Math.floor((double)iXYtl[1])/this.overlapStep))="+
(((Math.floor((double)iXYbr[1])/this.overlapStep))>((Math.floor((double)iXYtl[1])/this.overlapStep)))+
" "+((int) (Math.floor((double)iXYbr[1])/this.overlapStep))+">"+((int) (Math.floor((double)iXYtl[1])/this.overlapStep))+"="+
(((int) (Math.floor((double)iXYbr[1])/this.overlapStep))>((int) (Math.floor((double)iXYtl[1])/this.overlapStep)))*/
);
}
/*
((int) (Math.floor((double)iXYbr[0])/this.overlapStep))>((int) (Math.floor((double)iXYtl[0])/this.overlapStep)),
((int) (Math.floor((double)iXYbr[1])/this.overlapStep))>((int) (Math.floor((double)iXYtl[1])/this.overlapStep))};
*/
extraBits[index++]=enMasks[whereY[y][1]][whereX[x][1]][whereY[y][0]*this.paddedSize + whereX[x][0]];
}
}
if (debugLevel>5){
(new ShowDoubleFloatArrays()).showArrays(
extraBits,
longPaddedSize,
longPaddedSize,
"EB_-X"+tileX+"-Y"+tileY
);
}
int index=0;
double [] result=new double[this.paddedSize*this.paddedSize];
double [][]k={
{(1-dXY[0])*(1-dXY[1]), ( dXY[0])*(1-dXY[1])},
{(1-dXY[0])*( dXY[1]), ( dXY[0])*( dXY[1])}};
for (int y=0;y=this.tilesX) || (tileY>=this.tilesY) || (thisZMap[tileY][tileX]==null)){
boolean [] empty=new boolean[this.paddedSize*this.paddedSize];
for (int i=0;i5){
boolean [] geno=thisZMap[tileY][tileX].getEnabledNonOccluded(
disparityThresholdFg,
disparity,
disparityTolerance,
debugLevel);
int dbgSize=(int)Math.sqrt(geno.length);
(new ShowDoubleFloatArrays()).showArrays(
geno,
dbgSize,
dbgSize,
"GENO_-X"+tileX+"-Y"+tileY
);
}
return thisZMap[tileY][tileX].getEnabledNonOccluded(
disparityThresholdFg,
disparity,
disparityTolerance,
debugLevel);
}
/**
* Calculate fine-step (4 pix?) correlation between the pair of square images (32x32) that are already shifted to compensate for known disparity
* only center area (with 1/4 margins) will be calculated, if the correlation is strong enough, the updated disparity arrays will be calculated fro these tiles
*
* @param data a pair of square arrays to be correlated (32x32 pixels)
* @param window cosine mask (8x8=64 pixels long). If the center element is 0.0 - will generate array
* @param dXY unity vector in the direction of the disparity between the two images
* @param phaseCoeff 0.0 - normal correlation, 1.0 - pure phase correlation
* @param debugLevel debug level
* @return smooth array, same dimension as data[0], with only center area trusted, proportional to correlation at zero
*/
public double [] localCorrelation(
DoubleFHT doubleFHT, // will generate if null, can be used to reuse initialization
double [][] data, // difines size
double [] window, // defines tile size
double [] dXY, // unity vector defines disparity direction
double phaseCoeff, // phase correlation fraction (0 - normal, 1.0 - pure phase)
double highPassSigma,
double lowPassSigma,
int normalize,
int debugLevel
){
return localCorrelation(
doubleFHT, // will generate if null, can be used to reuse initialization
data, // difines size
window, // defines tile size
dXY, // unity vector defines disparity direction
phaseCoeff, // phase correlation fraction (0 - normal, 1.0 - pure phase)
highPassSigma,
lowPassSigma,
normalize,
debugLevel,
null);
}
public double [] localCorrelation(
DoubleFHT doubleFHT, // will generate if null, can be used to reuse initialization
double [][] data, // difines size
double [] window, // defines tile size
double [] dXY, // unity vector defines disparity direction
double phaseCoeff, // phase correlation fraction (0 - normal, 1.0 - pure phase)
double highPassSigma,
double lowPassSigma,
int normalize,
int debugLevel,
String title
){
int size=(int) Math.sqrt(data[0].length); // 32
int tileSize=(int) Math.sqrt(window.length); // 8
int tileStep=tileSize/2; //4
int tileExtra=size-tileSize; // 24??
int tileLength=tileSize*tileSize; // 64
int tileCenter=(tileSize+1)*tileSize/2; // 36
int numTilesRow=size/(2*tileStep); // 4
if (doubleFHT==null) doubleFHT= new DoubleFHT();
if (Double.isNaN(window[0])){
double [] window1d=doubleFHT.getHamming1d(tileSize); // Hamming
for (int i=0;i3) && (tileY==0) && (tileX==0)) for (int i=0;i4){
double d=Math.sqrt(dxA*dxA+dyA*dyA);
System.out.println("getShiftedSlices(): xc1="+xc1+" yc1="+yc1+" size="+size+" first="+first+" second="+second);
System.out.println("getShiftedSlices(): d="+IJ.d2s(d,2)+" dX="+IJ.d2s(dxA,2)+" dY="+IJ.d2s(dyA,2));
System.out.println("getShiftedSlices(): shiftImage[0]="+shiftImage[0]+" shiftImage[1]="+shiftImage[1]);
System.out.println("getShiftedSlices() first image: dX="+IJ.d2s(dx[0],2)+" dY="+IJ.d2s(dy[0],2));
System.out.println("getShiftedSlices() second image: dX="+IJ.d2s(dx[1],2)+" dY="+IJ.d2s(dy[1],2));
System.out.println("getShiftedSlices() second image: iDx="+iDx+" iDyY="+iDy);
System.out.println("getShiftedSlices() second image: xc="+(xc1+iDx)+" yc="+(yc1+iDy));
}
int doubleSize=2*size;
int doubleLength=doubleSize*doubleSize;
double [] slice= new double[doubleLength];
double [] imgSlice;
int [] img= {first,second};
/*
if (window==null) {
window = new double [doubleLength];
window[0]=Double.NaN;
}
*/
if (Double.isNaN(window[0])) {
double[] window1d=doubleFHT.getHamming1d(size); // Hamming
int index=0;
int halfSize=size/2;
int size32=3*size/2;
for (int iy=0;iysize32)?window1d[iy-size]:1.0);
for (int ix=0;ixsize32)?window1d[ix-size]:1.0);
window[index++]=wx*wy;
}
}
}
if (debugLevel>4){
System.out.print("getShiftedSlices(): channelMask="+channelMask);
if (debugLevel>3){
(new ShowDoubleFloatArrays()).showArrays(
window,
doubleSize,
doubleSize,
"window");
}
}
for (int iImg=0;iImg2) System.out.println("xc["+iImg+"]="+xc[iImg]+" yc["+iImg+"]="+yc[iImg]+" shiftImage["+iImg+"]="+shiftImage[iImg]);
slice= getSelection(
imgSlice, // source image/channel slice
slice,
doubleSize, //int width,
doubleSize, //int height,
imageWidth,
xc[iImg] , //xc,
yc[iImg]); //yc);
if (debugLevel>4) {
(new ShowDoubleFloatArrays()).showArrays(
slice,
doubleSize,
doubleSize,
"slice-getShiftedSlices");
}
// normalize and save DC
double dc=0.0;
if (shiftImage[iImg]){
dc=normalizeAndWindowGetDC (slice, window); //windowInterpolation
if (debugLevel>4) {
(new ShowDoubleFloatArrays()).showArrays(
slice,
doubleSize,
doubleSize,
"slice-normalized-dc"+dc);
}
// doubleFHT.shift(slice, dx[iImg], dy[iImg]);
doubleFHT.shift(slice, -dx[iImg], -dy[iImg]);
if (debugLevel>4) {
(new ShowDoubleFloatArrays()).showArrays(
slice,
doubleSize,
doubleSize,
"slice-shifted");
}
} else {
if (debugLevel>2) System.out.println("No shift is needed");
}
// int oSliceNumber=img.length*cN+iImg;
if (debugLevel>2) System.out.println("getShiftedSlices() dc="+dc);
if (doubleSizeOutput) {
if (debugLevel>2) System.out.println("doubleLength="+doubleLength+"doubleSize="+doubleSize);
result[channel][iImg]=new double [doubleLength];
int oIndex=0;
for (int iY=0;iY=result[channel][iImg].length) || (oIndex>=slice.length) || (oIndex>=window.length)){
System.out.println("iX="+iX+" iY="+iY+" oIndex="+oIndex);
}
result[channel][iImg][oIndex]=shiftImage[iImg]?(slice[oIndex]/window[oIndex]+dc):slice[oIndex];// no NaN with Hamming
oIndex++;
}
}
}else{
result[channel][iImg]=new double [size*size];
int iIndex=size*size+size/2; // top left index of the inner square size*size
int oIndex=0;
for (int iY=0;iY=fullHeight);
for (int ix=0;ix=fullWidth)) {
if (oIndex>=selection.length) System.out.println("\ngetSelection(imageSlice["+imageSlice.length+"],selection["+selection.length+"]"+
","+fullWidth+","+width+","+height+","+xc+","+yc+") ix="+ix+" iy="+iy+" oIndex="+oIndex);
selection[oIndex]=0.0;
} else {
selection[oIndex]=imageSlice[srcY*fullWidth+srcX];
}
}
}
return selection;
}
public void setupZMap(
Rectangle woi, // in tiles - may be
final int nImg,
final int maxNumber,
final double minFirst,
final double minAbsolute,
final double minRelative,
final double mergeMax,
final int overlap,
final double zMapMinForeground,
final int threadsMax,
final boolean showProgress,
final int debugLevel){
if (this.centerPixels==null) {
String msg="Centered disparity data is not defined";
IJ.showMessage("Error",msg);
throw new IllegalArgumentException (msg);
}
if (this.syntheticPixels==null) moveDisparityFromCenter(threadsMax,showProgress,debugLevel);
final float [] thisSyntheticPixels= this.syntheticPixels[nImg];
if (woi==null) woi=new Rectangle (0, 0, this.tilesX,this.tilesY);
this.zMapWOI=new Rectangle();
this.zMapWOI.x=(woi.x>=0)?woi.x:0;
this.zMapWOI.y=(woi.y>=0)?woi.y:0;
this.zMapWOI.width= (((woi.x+woi.width) <=this.tilesX)?(woi.x+woi.width): this.tilesX)-this.zMapWOI.x;
this.zMapWOI.height=(((woi.y+woi.height)<=this.tilesY)?(woi.y+woi.height):this.tilesY)-this.zMapWOI.y;
final Rectangle zMapWOI=this.zMapWOI;
if ((this.zMap==null) || (this.zMap[nImg]==null)) initZMap(nImg);
final ZTile [][] thisZMap=this.zMap[nImg];
final int tilesX=this.zMapWOI.width;
final int tiles=this.zMapWOI.width*this.zMapWOI.height;
if (debugLevel>2) System.out.println("setupZMap() woi.x="+woi.x+" woi.y="+woi.y+
" woi.width="+woi.width+" woi.height="+woi.height);
if (debugLevel>2) System.out.println("setupZMap() zMapWOI.x="+zMapWOI.x+" zMapWOI.y="+zMapWOI.y+
" zMapWOI.width="+zMapWOI.width+" zMapWOI.height="+zMapWOI.height);
this.paddedSize=this.overlapStep+2*overlap;
this.innerMask=new BitSet(this.paddedSize*this.paddedSize);
for (int i=0;ithis.tilesY) tileYmax=this.tilesY;
int tileXmin=woi.x/this.overlapStep;
int tileXmax=(xRight-1)/this.overlapStep+1; // +1 to the last needed
if (tileXmax>this.tilesX) tileXmax=this.tilesX;
for (int i=0;i=0) && (oY=0) && (oXthis.tilesY) tileYmax=this.tilesY;
int tileXmin=woi.x/this.overlapStep;
int tileXmax=(xRight-1)/this.overlapStep+1; // +1 to the last needed
if (tileXmax>this.tilesX) tileXmax=this.tilesX;
for (int i=0;i=0) && (oY=0) && (oX2) System.out.println("setupZMapTile("+tileX+","+tileY+",...)");
boolean debugThisTile=debugLevel>3;
thisZMap[tileY][tileX]=new ZTile();
ZTile zTile=thisZMap[tileY][tileX];
double [] correlationSection = getTileDisparity(tileX, tileY, thisDisparityTiles);
// double threshold=minAbsolute;
boolean [] isMax =new boolean[correlationSection.length];
int iMax=-1;
int numMax=0;
for (int i=0;i=minAbsolute) &&
((i==0)||(correlationSection[i]>=correlationSection[i-1])) &&
((i==(correlationSection.length-1)) || (correlationSection[i]>=correlationSection[i+1]))){
isMax[i]=true; // local max above absolute threshold
if ((numMax==0) || (correlationSection[i]>correlationSection[iMax])) iMax=i;
numMax=1;
} else {
isMax[i]=false;
}
}
// double mergeMax=1.5; //merge maximums for interpolation if closer than mergeMax;
int iMergeMax=(int) Math.round(mergeMax/this.disparityPerEntry);
if (debugThisTile){
System.out.println("mergeMax="+mergeMax+" iMergeMax="+iMergeMax);
}
if ((numMax>0) && (correlationSection[iMax]=threshold) numMax++;
numMax++;
//else isMax[i]=false;
}
}
if (numMax>maxNumber) numMax=maxNumber; // limit to specified number of correlation maximums to subpixel
int [] maxIndices=new int [numMax];
maxIndices[0]=iMax;
isMax[iMax]=false;
int maxNum;
for (maxNum=1;maxNum=isMax.length) break;
if (isMax[maxIndices[maxNum-1]+i]){
nearUp = maxIndices[maxNum-1]+i;
break;
}
}
int n=1+((nearDown>=0)?1:0)+((nearUp>=0)?1:0);
if (n>1){
double s=maxIndices[maxNum-1]+((nearDown>=0)?nearDown:0)+((nearUp>=0)?nearUp:0);
int iMerged=(int) Math.round(s/n);
if (debugThisTile){
System.out.println("Merging close maximums: "+maxIndices[maxNum-1]+" with "+
((nearDown>=0)?nearDown:"")+" "+((nearUp>=0)?nearUp:"")+" to "+iMerged);
}
maxIndices[maxNum-1]=iMerged;
isMax[iMerged]=false;
if (nearDown>=0) isMax[nearDown]=false;
if (nearUp>=0) isMax[nearUp]=false;
}
iMax=-1;
for (int i=0;i=correlationSection[iMax]))) iMax=i;
}
if (iMax<0) break; // no more maximums
isMax[iMax]=false;
maxIndices[maxNum]=iMax;
}
if (debugThisTile){
System.out.println("List maximums on correlation section");
for (int n=0;n=this.enabledPixels.length) return; // or make it an error?
int paddedSize= this.size+this.overlap*2;
this.enabledPixels[this.foregroundIndex]=new BitSet(paddedSize*paddedSize);
this.enabledPixels[this.foregroundIndex].set(0,paddedSize*paddedSize);
this.enabledPlane[this.foregroundIndex]=true;
}
public boolean isForegroundValid(){
if (this.foregroundIndex>=this.maximums.length) return false; // no foreground
if (!this.enabledPlane[this.foregroundIndex]) return false; // disabled plane
if (this.enabledPixels[this.foregroundIndex]== null) return true; // all pixels enabled
BitSet bb=(BitSet) this.enabledPixels[this.foregroundIndex].clone();
bb.and(this.innerMask);
if (!bb.isEmpty()) return true;
return false;
}
public boolean advanceForeground(){
while (!isForegroundValid()){
if (this.foregroundIndex>=this.maximums.length) return false; // no foreground left
this.foregroundIndex++;
}
return true;
}
public float [] getZmap(){
if (this.zmap!=null) return this.zmap;
int paddedSize=this.size+2*this.overlap;
int paddedLength=paddedSize*paddedSize;
float [] map= new float [paddedLength];
BitSet needed=null;
BitSet newBits=null;
for (int plane=this.foregroundIndex; plane=0;i=newBits.nextSetBit(i+1)){
map[i]=(float) this.maximums[plane][0];
}
needed.andNot(newBits);
if (this.enabledPixels[plane]==null) return map;
if (needed.isEmpty()) return map;
}
if (needed==null){ // no maximums at all, probably
for (int i=0;i=0;i=needed.nextSetBit(i+1))map[i]=0.0F;
return map;
}
/*
public boolean [] getNonOccluded(double disparity){
int paddedSize=this.size+2*this.overlap;
int paddedLength=paddedSize*paddedSize;
BitSet nonOccludedBits=new BitSet(paddedLength);
nonOccludedBits.set(0,paddedLength);
for (int plane=0;(planedisparity);plane++) if (this.enabledPlane[plane]) {
if (this.certainPixels[plane]!=null) nonOccludedBits.andNot(this.certainPixels[plane]);
}
boolean [] nonOccluded=new boolean[paddedLength];
for (int i=0;i3) System.out.println(" ---- getEnabledNonOccluded("+disparityFG+","+
disparity+","+
disparityTolerance+","+
debugLevel+")");
int paddedSize=this.size+2*this.overlap;
int paddedLength=paddedSize*paddedSize;
BitSet nonOccludedBits=new BitSet(paddedLength);
nonOccludedBits.set(0,paddedLength);
for (int plane=0;(planedisparityFG);plane++) if (this.enabledPlane[plane]) {
if (this.certainPixels[plane]!=null) nonOccludedBits.andNot(this.certainPixels[plane]);
if (debugLevel>3) System.out.println(" ------ plane="+plane+" cumul. nonOccludedBits.cardinality()="+nonOccludedBits.cardinality());
}
if (!Double.isNaN(disparityTolerance)){
BitSet enabledBits=new BitSet(paddedLength);
for (int plane=0;plane3) System.out.println(" ------ plane="+plane+" cumul. enabledBits.cardinality()="+enabledBits.cardinality());
}
nonOccludedBits.and(enabledBits);
if (debugLevel>3) System.out.println(" ------ result nonOccludedBitscardinality()="+nonOccludedBits.cardinality());
}
boolean [] nonOccluded=new boolean[paddedLength];
for (int i=0;iaMax) aMax=this.maximums[i][1];
aMax=Math.max(aMax*this.minRelative,this.minAbsolute);
for (int i=0;i=aMax;
if (keepFG && (this.foregroundIndex=0)) {
this.enabledPlane[this.foregroundIndex]=true;
}
}
public void setForeGround(
double minForeground){
if (!Double.isNaN(minForeground)) this.foregroundThreshold=minForeground;
for (this.foregroundIndex=0;
(this.foregroundIndex=this.subdivAverage) countAvg=this.subdivAverage-1;
int countDiff= (int)Math.round(0.5*(v2-v1))+subdivHalfDifference;
if (countDiff<0) countDiff=0;
if (countDiff>2*this.subdivHalfDifference) countDiff=2*this.subdivHalfDifference;
staging[countDiff*this.subdivAverage+countAvg]+=weight;
}
public void addStagingSample(
double []staging,
int chn,
double weight,
double v1,
double v2,
// reduce weight depending on difference, scale to measured variance
double scaleVariance,
double kLocal, // 0 - use global varaiance, 1.0 - use local
int nImg1, // first image number
int nImg2, // second image number
double var1, // first image variance
double var2, // second image variance
boolean debug){
double min=this.valueLimits[chn][0];
double step=(this.valueLimits[chn][1]-this.valueLimits[chn][0])/(this.subdivAverage-1);
double av=0.5*(v1+v2);
int countAvg=(int) Math.round((av-min)/step);
if (countAvg<0) countAvg=0;
if (countAvg>=this.subdivAverage) countAvg=this.subdivAverage-1;
double diff=0.5*(v2-v1);
int countDiff= (int)Math.round(diff/step)+subdivHalfDifference;
if (countDiff<0) countDiff=0;
if (countDiff>2*this.subdivHalfDifference) countDiff=2*this.subdivHalfDifference;
if (debug){
System.out.println("addStagingSample(...,"+chn+","+weight+","+v1+","+v2+","+scaleVariance+","+
kLocal+","+nImg1+","+nImg2+","+var1+","+var2);
System.out.println(" +++ min="+min+" step="+step+" av="+av+" diff="+diff+" countAvg="+countAvg+" countDiff="+countDiff);
}
if (scaleVariance>0.0){
if (kLocal<1.0){ // mix local variance with average over all image
int count1=(int) Math.round((v1-min)/step);
if (count1<0) count1=0;
if (count1>=this.subdivAverage) count1=this.subdivAverage-1;
int count2=(int) Math.round((v2-min)/step);
if (count2<0) count2=0;
if (count2>=this.subdivAverage) count2=this.subdivAverage-1;
if (kLocal<0) kLocal=0;
var1=kLocal*var1+(1.0-kLocal)*this.standardVariance[nImg1][chn][count1];
var2=kLocal*var2+(1.0-kLocal)*this.standardVariance[nImg2][chn][count2];
}
double sigma=scaleVariance*Math.sqrt(var1*var1);
weight*=Math.exp(-diff*diff/(2.0*sigma*sigma))/sigma/Math.sqrt(2*Math.PI);
if (debug) System.out.println(" +++ sigma="+sigma+" weight="+weight);
}
staging[countDiff*this.subdivAverage+countAvg]+=weight;
if (debug) System.out.println(" staging["+(countDiff*this.subdivAverage+countAvg)+"]="+staging[countDiff*this.subdivAverage+countAvg]);
}
public double getStrength(
double []staging,
int chn,
double v1,
double v2){
double min=this.valueLimits[chn][0];
double step=(this.valueLimits[chn][1]-this.valueLimits[chn][0])/(this.subdivAverage-1);
double av=0.5*(v1+v2);
int countAvg=(int) Math.floor((av-min)/step);
int countAvg1=(int) Math.ceil((av-min)/step);
double dx=(av-min)-step*countAvg;
if (countAvg<0) {
countAvg=0;
countAvg1=0;
dx=0.0;
}
if (countAvg1>=this.subdivAverage){
countAvg=this.subdivAverage-1;
countAvg1=this.subdivAverage-1;
dx=0.0;
}
double diff=0.5*(v2-v1);
int countDiff= (int)Math.floor(diff/step)+this.subdivHalfDifference;
int countDiff1= (int)Math.ceil(diff/step)+this.subdivHalfDifference;
double dy=diff-step*(countDiff-this.subdivHalfDifference);
if (countDiff<0) {
countDiff=0;
countDiff1=0;
dy=0.0;
}
if (countDiff1>2*this.subdivHalfDifference){
countDiff=2*this.subdivHalfDifference;
countDiff1=2*this.subdivHalfDifference;
dy=0.0;
}
if ((countDiff1*this.subdivAverage+countAvg1)>=staging.length){
System.out.println("BUG: getStrength() countAvg="+countAvg+" countAvg1="+countAvg1+
" countDiff="+countDiff+" countDiff1="+countDiff1+" staging.length="+staging.length);
}
return
staging[countDiff *this.subdivAverage+countAvg ]*(1-dy)*(1-dx)+
staging[countDiff *this.subdivAverage+countAvg1]*(1-dy)*( dx)+
staging[countDiff1*this.subdivAverage+countAvg ]*( dy)*(1-dx)+
staging[countDiff1*this.subdivAverage+countAvg1]*( dy)*( dx);
}
public void blurStaging(
double []staging,
int chn,
double sigma){
if (sigma<=0) return;
// double min=this.valueLimits[chn][0];
double step=(this.valueLimits[chn][1]-this.valueLimits[chn][0])/(this.subdivAverage-1);
(new DoubleGaussianBlur()).blurDouble(
staging,
this.subdivAverage,
2*this.subdivHalfDifference+1,
sigma/step,
sigma/step,
0.01);
}
public Photometric(
double [][][] images,
int imageWidth,
int margins,
double ignoreFraction,
int subdivAverage,
int subdivHalfDifference,
// int subdivVariance,
double smoothVarianceSigma,
double scaleVariance,
int debugLevel
){
if (margins<1) margins=1;
this.subdivAverage=subdivAverage;
this.subdivHalfDifference=subdivHalfDifference;
// this.subdivVariance = subdivVariance;
this.smoothVarianceSigma=smoothVarianceSigma;
this.scaleVariance=scaleVariance; // sigma along difference (vertical axis) will be scaleVariance*variance for this average value
this.numImages=images.length;
this.valueLimits=new double[channelNames.length][];
int [] dirs1={1,imageWidth+1,imageWidth,imageWidth-1,-1,-imageWidth-1,-imageWidth,-imageWidth+1,0};
this.standardVariance=new double[this.numImages][this.channelNames.length][];
this.averageVariance=new double[this.numImages][this.channelNames.length];
for (int i=0;i2){
System.out.println("nImg="+nImg+" chn="+chn+
" min0="+IJ.d2s(min,3)+" max0="+IJ.d2s(max,3));
}
for (int y=margins;y<(imageHeight-margins);y++) for (int x=margins;x<(imageWidth-margins);x++){
if (image[y*imageWidth+x]2) System.out.println("y="+y+" x="+x+" min="+IJ.d2s(min,3));
}
if (image[y*imageWidth+x]>max) {
max=image[y*imageWidth+x];
if (debugLevel>2) System.out.println("y="+y+" x="+x+" max="+IJ.d2s(max,3));
}
}
if (debugLevel>2){
System.out.println("nImg="+nImg+" chn="+chn+
" min="+IJ.d2s(min,3)+" max="+IJ.d2s(max,3));
}
}
int [] histogram=new int [this.histogramSize];
for (int i=0;i2){
System.out.println( " min="+IJ.d2s(min,3)+" max="+IJ.d2s(max,3)+" step="+IJ.d2s(step,6));
}
for (int nImg=0;nImg=this.histogramSize) index=this.histogramSize-1;
histogram[index]++; //java.lang.ArrayIndexOutOfBoundsException: 1005
}
}
int totalNum=0;
for (int i=0;i2){
System.out.println( " totalNum="+totalNum+" ignoreNumPix="+ignoreNumPix);
}
this.valueLimits[chn]=new double[2];
int num=0;
for (int i=0;i2) System.out.println("--- "+num+" min+"+i+"*step="+IJ.d2s(min+i*step,3));
if (num>=ignoreNumPix){
this.valueLimits[chn][0]=min+i*step;
if (debugLevel>2){
System.out.println("i="+i+" min+i*step="+IJ.d2s(min+i*step,3));
}
break;
}
}
num=0;
for (int i=histogram.length-1;i>=0;i--){
num+=histogram[i];
if (debugLevel>2) System.out.println("--- "+num+" min+"+i+"*step="+IJ.d2s(min+i*step,3));
if (num>=ignoreNumPix){
this.valueLimits[chn][1]=min+i*step;
if (debugLevel>2){
System.out.println("i="+i+" min+i*step="+IJ.d2s(min+i*step,3));
}
break;
}
}
if (debugLevel>1){
System.out.println("Channel '"+this.channelNames[chn]+
"' min="+IJ.d2s(this.valueLimits[chn][0],3)+" max="+IJ.d2s(this.valueLimits[chn][1],3));
}
// calculatye variance for each emage/channel
for (int nImg=0;nImg=this.subdivAverage) count=this.subdivAverage-1;
this.standardVariance[nImg][chn][count]+=v2; // squared
samples[count]++;
this.averageVariance[nImg][chn]+=v2;
totalSamples++;
}
for (int i=0;i0) this.standardVariance[nImg][chn][i]=Math.sqrt(this.standardVariance[nImg][chn][i]/samples[i]);
else this.standardVariance[nImg][chn][i]=0.0;
// if ((debugLevel>1) && (nImg==0) && (chn==0)){
if ((debugLevel>2) && (nImg==0)){
System.out.println("samples["+i+"]="+samples[i]+
" this.standardVariance["+nImg+"]["+chn+"]["+i+"]="+IJ.d2s(this.standardVariance[nImg][chn][i],3));
}
}
this.averageVariance[nImg][chn]=Math.sqrt(this.averageVariance[nImg][chn]/totalSamples);
for (int i=0;i0.0){
(new DoubleGaussianBlur()).blur1Direction(
this.standardVariance[nImg][chn], //double [] pixels,
this.subdivAverage, //int width,
1, //int height,
this.smoothVarianceSigma, //double sigma,
0.01, //double accuracy,
true //boolean xDirection: true - horizontal
);
}
}
}
// now cteate matchingQuality arrays for each image pair/channel
// this.matchingQuality=new double [this.numImages][this.numImages-1][this.valueLimits.length][this.subdivAverage*this.subdivDifference];
this.matchingQuality=new double [this.numImages][this.numImages-1][][]; //[this.subdivAverage*this.subdivDifference];
int subdivDifference=2*this.subdivHalfDifference+1;
int matchingQualityLength=this.subdivAverage*subdivDifference;
for (int nImg=0;nImg=nImg)?(sIndex+1):sIndex;
if ((images[nImg]!=null) && (images[sImg]!=null)){
this.matchingQuality[nImg][sIndex]=new double [this.valueLimits.length][];
for (int chn=0;chn=nImg)?(sIndex+1):sIndex;
for (int chn=0;chn2) debugArray=new double[3][data[0].length*subPixel*subPixel];
for (int tileY=0;tileY3) && (tileY==0) && (tileX==0)) for (int i=0;i2){
System.out.println("calculateDisparityTiles(), tileXMin="+tileXMin+", tileXMax="+tileXMax+", tileYMin="+tileYMin+", tileYMax="+tileYMax);
}
final double [][][][][] disparityTiles= new double [numTilesY][numTilesX][][][]; // first [tileY][tileX][image number][other image index][disparity index]
final double [][][][] fhtCache=new double [numTilesY][][][];
for (int tileY=tileYMin;tileY<=tileYMax;tileY++){
for (int tileX=tileXMin;tileX<=tileXMax;tileX++){
disparityTiles[tileY][tileX]=null;
}
fhtCache[tileY]=null;
}
final int numSensors=this.channel.length;
int minDispalcementTileXAll=minDispalcementTileX[0];
int maxDispalcementTileXAll=maxDispalcementTileX[0];
int minDispalcementTileYAll=minDispalcementTileY[0];
int maxDispalcementTileYAll=maxDispalcementTileY[0];
for (int i=1;iminDispalcementTileX[i]) minDispalcementTileXAll=minDispalcementTileX[i];
if (maxDispalcementTileXAllminDispalcementTileY[i]) minDispalcementTileYAll=minDispalcementTileY[i];
if (maxDispalcementTileYAll2){
System.out.println("calculateDisparityTiles(), minDispalcementTileXAll="+minDispalcementTileXAll+", maxDispalcementTileXAll="+maxDispalcementTileXAll);
System.out.println("calculateDisparityTiles(), minDispalcementTileXAll="+minDispalcementTileXAll+", maxDispalcementTileXAll="+maxDispalcementTileXAll);
System.out.println("calculateDisparityTiles(), minDispalcementTileYAll="+minDispalcementTileYAll+", maxDispalcementTileYAll="+maxDispalcementTileYAll);
System.out.println("calculateDisparityTiles(), minDispalcementTileYAll="+minDispalcementTileYAll+", maxDispalcementTileYAll="+maxDispalcementTileYAll);
for (int i=0;i2){
System.out.println("==== calculateDisparityTiles(), correlating channel "+channelName[finalChnYCbCr]+" tileY="+tileY+" ====");
}
// free memory - delete the full row of unneeded FHT tiles (all images)
if ((tileY+minDispalcementTileYAll)>0) {
if (debugLevel>2){
System.out.println("calculateDisparityTiles(), freeing tile row "+(tileY+minDispalcementTileYAll-1)+" for all images");
}
fhtCache[tileY+minDispalcementTileYAll-1]=null;
for (int nImg=0;nImg=0){
if (debugLevel>2){
System.out.println("calculateDisparityTiles(), freeing tile row "+tileYDel+" for image "+nImg);
}
if (fhtCache[tileYDel]!=null) fhtCache[tileYDel][nImg]=null;
if (nextCachedY[nImg]<=tileYDel) nextCachedY[nImg]= tileYDel+1; // for filling up cache initially only
}
}
// calculate rows of FHT tiles
for (int nImg=0;nImg=numTilesX)correctedTileXMax=numTilesX-1;
for (;nextCachedY[nImg]<=tileY+maxDispalcementTileY[nImg];nextCachedY[nImg]++) if (nextCachedY[nImg]=0.0) || (correlationLowPassSigma>=0.0);
if (useFilter) doubleFHT.createFrequencyFilter(
window, // just to specify size
correlationHighPassSigma,
correlationLowPassSigma);
double [][][] stage=null; // [pY][pX]{data, weight}
for (int rsltTile=rsltTileAtomic.getAndIncrement(); rsltTile2) && (finalChnYCbCr==0) && (tileY==tileYMin) &&(rsltTile==0); // make
boolean debugThisTile=(debugLevel>2) && (finalChnYCbCr==0) && (tileY==tileYMin) &&(rsltTile==1); // make
int nImg=firstSecondIndex[rsltTile%firstSecondIndex.length][0];
int secondIndex=firstSecondIndex[rsltTile%firstSecondIndex.length][1];
int secondImage=otherImage[nImg][secondIndex];
int [][] thisDisplacementTilesList=displacementTilesList[nImg][secondIndex];
if ((stage==null) ||
(stage.length!= correlationMap[nImg][secondIndex].length) ||
(stage[0].length!=correlationMap[nImg][secondIndex][0].length) ){
stage=new double [correlationMap[nImg][secondIndex].length][correlationMap[nImg][secondIndex][0].length][2];
}
double [] zero2={0.0,0.0};
// TODO - zero only needed pixels?
for (int pY=0;pY=0) && (tileY=0) && (tileX=0) && (tY