codebase / stimlib

Commit a47a23a9dcd413b69032b0c4a5b74c72db664fb8

Authored by dmayerich
1 parent f8303137

added ENVI functions

Showing 17 changed files with 1825 additions and 208 deletions   Show diff stats
Inline Side-by-side
envi/EnviClose.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef ENVICLOSE_H
  2 +#define ENVICLOSE_H
  3 +
  4 +
  5 +#include <fstream>
  6 +#include <stdio.h>
  7 +#include "EnviHeader.h"
  8 +#include "EnviFid.h"
  9 +#include "EnviOpen.h"
  10 +
  11 +// this is to close all the opened files and to create a header file for the written files
  12 +
  13 +static void EnviClose(EnviFidClass &fid, int dimX=0)
  14 +{
  15 +
  16 + float* zeroPtr = (float*) malloc(sizeof(float) * fid.header.nF);
  17 + memset(zeroPtr, 0, sizeof(float) * fid.header.nF);
  18 +
  19 + if(fid.mode=="r")
  20 + {
  21 + fclose(fid.file);
  22 + }
  23 +
  24 + else
  25 + {
  26 + int totalPixels = fid.header.sX;
  27 + int dimY = totalPixels/dimX;
  28 +
  29 + if(totalPixels % dimX > 0)
  30 + {
  31 + dimY++;
  32 +
  33 + int padding = dimY*dimX-totalPixels;
  34 +
  35 + for(int p=0; p < padding; p++)
  36 + fwrite(zeroPtr, sizeof(float), fid.header.nF, fid.file);
  37 + }
  38 +
  39 + //write the calculated values to the header object
  40 + fid.header.sX = dimX;
  41 + fid.header.sY = dimY;
  42 +
  43 +
  44 + //write the header file to disk
  45 + //TODO
  46 +
  47 +
  48 + //close the binary file
  49 + fclose(fid.file);
  50 +
  51 + }
  52 +
  53 +
  54 +
  55 +
  56 +}
  57 +#endif
  58 +
... ...
envi/EnviFid.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef ENVIFID_H
  2 +#define ENVIFID_H
  3 +
  4 +#include <string>
  5 +#include <iostream>
  6 +#include <fstream>
  7 +#include <ctime>
  8 +#include <sstream>
  9 +#include "EnviHeader.h"
  10 +using namespace std;
  11 +
  12 +class EnviFidClass
  13 +{
  14 +public:
  15 + //binary file name
  16 + string fileName;
  17 +
  18 + //header structure
  19 + EnviHeaderClass header;
  20 +
  21 + //valid flag (is the current file valid)
  22 + int isValid;
  23 +
  24 + //file mask (NULL if the whole file is to be read)
  25 + char* mask;
  26 +
  27 + //current position (in pixels) in the file
  28 + int filePos;
  29 +
  30 + //read/write mode ("r" or "w")
  31 + string mode;
  32 +
  33 + //create a variable to store a pointer to the binary file
  34 + FILE * file;
  35 +
  36 + //constructor, set default values
  37 + EnviFidClass()
  38 + {
  39 + isValid = 0;
  40 + filePos = 0;
  41 + }
  42 +
  43 +
  44 +};
  45 +
  46 +
  47 +#endif
... ...
envi/EnviHeader.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef ENVIHEADER_H
  2 +#define ENVIHEADER_H
  3 +#include <string>
  4 +#include <iostream>
  5 +#include <fstream>
  6 +#include <ctime>
  7 +#include <sstream>
  8 +using namespace std;
  9 +
  10 +enum InterleaveType {bip, bil, bsq};
  11 +class EnviHeaderClass //"EnviheaderClass" will be used to store data into variables sX, sY, etc.
  12 +{
  13 +public:
  14 + unsigned int sX; //samples = "sX"
  15 + unsigned int sY; //lines = "sY"
  16 + unsigned int nF; //bands = "nF"
  17 + unsigned int headersize; //header offset = "headersize"
  18 + unsigned int datatype; //data type = "datatype"
  19 + InterleaveType type;
  20 + string units; //wavelength units = "units" (Usully in Wavenumber or Wavelength)
  21 + int LoadHeader(string filename)
  22 + {
  23 + ifstream inFile;
  24 + //Load the header file in the same directory, for future note, might want to modify code to either "cin" the file name of have a feature in the GUI which take into account different fileneames
  25 + inFile.open (filename.c_str());
  26 +
  27 + //If or not there is a file
  28 + if (!inFile)
  29 + {
  30 + cout<<"Error"<<endl;
  31 + return 1;
  32 + }
  33 +
  34 + string line;
  35 + getline(inFile,line);
  36 +
  37 + //If the 1st line of text file does not read ENVI, spit out an error.
  38 + if (line.compare("ENVI"))
  39 + {
  40 + cout<<"Error, not an ENVI file"<<endl;
  41 + return 1;
  42 + }
  43 +
  44 + while(inFile) //Loops through entire file
  45 + {
  46 + getline(inFile,line);
  47 + if(line.find("samples")==0)
  48 + {
  49 + sX=getNumber(line); //Calls getNumber function for "samples", stores in header.sX
  50 + }
  51 + if(line.find("lines")==0)
  52 + {
  53 + sY=getNumber(line); //Calls getNumber function for "lines", stores in header.sY
  54 + }
  55 + if(line.find("bands")==0)
  56 + {
  57 + nF=getNumber(line); //Calls getNumber function for "bands", stores in header.nF
  58 + }
  59 + if(line.find("header offset")==0)
  60 + {
  61 + headersize=getNumber(line); //Calls getNumber function for "header offset", stores in header.headersize
  62 + }
  63 + if(line.find("data type")==0)
  64 + {
  65 + datatype=getNumber(line); //Calls getNumber function for "data type", stores in header.datatype
  66 + }
  67 +
  68 + line.find("interleave");
  69 + if(line.find("interleave")==0)
  70 + {
  71 + string typeStr;
  72 + typeStr=getString(line);
  73 + if (typeStr.find("bip")==0)
  74 + type=bip;
  75 + else if(typeStr.find("bsq")==0)
  76 + type=bsq;
  77 + else if(typeStr.find("bil")==0)
  78 + type=bil;
  79 + else
  80 + {
  81 + //throw an error and exit
  82 + cout<<"Type not recognized"<<endl;
  83 + return 1;
  84 + }
  85 +
  86 + }
  87 + if(line.find("wavelength units")==0)
  88 + {
  89 + //Calls getString function for "wavelength units", stores in header.units
  90 + units=getString(line);
  91 + }
  92 + }
  93 + return 0;
  94 + }
  95 +
  96 +
  97 +private:
  98 + //retrieves a numerical parameter from a line of text
  99 + int getNumber(string line)
  100 + {
  101 + int result=0;
  102 + int position;
  103 + position=line.find("="); //Finds the position of the equal sign
  104 +
  105 + string numStr; //Creates a string for storing data called "numStr"
  106 + numStr=line.substr(position+2); //Moves from the position of "=" 2 spaces to the right...
  107 + istringstream convert(numStr); //...and converts the string into numbers using "istringstream"
  108 +
  109 + convert>>result; //converts result from previous line
  110 +
  111 + return result;
  112 + }
  113 + //retrieves a string parameter from a line of text
  114 + string getString(string line)
  115 + {
  116 + int position;
  117 + position=line.find("="); //Finds the position of the equal sign
  118 +
  119 + string word; //Creates a string for storing word data called "word"
  120 + word=line.substr(position+2); //Moves from the position of "=" 2 spaces to the right
  121 + //Note: We do not need to convert word using "istringstream"
  122 +
  123 + return word;
  124 + }
  125 +};
  126 +#endif
... ...
envi/EnviOpen.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef ENVIOPEN_H
  2 +#define ENVIOPEN_H
  3 +#include <string>
  4 +#include <iostream>
  5 +#include <fstream>
  6 +#include <ctime>
  7 +#include <sstream>
  8 +#include <stdio.h>
  9 +#include "EnviHeader.h"
  10 +#include "EnviFid.h"
  11 +using namespace std;
  12 +
  13 +
  14 +
  15 +//return value is incorrect
  16 +static EnviFidClass EnviOpen(string filename, char* mask = NULL, string mode = "r", int nBands = 0, int sX = 0)
  17 +{
  18 + //create a file ID object
  19 + EnviFidClass fid;
  20 +
  21 + //set the file read/write mode
  22 + fid.mode = mode;
  23 +
  24 + //set the mask (if specified)
  25 + fid.mask = mask;
  26 +
  27 + // load the header file only while reading
  28 + if(mode=="r")
  29 + {
  30 + //load the header file
  31 + int error = fid.header.LoadHeader(filename);
  32 +
  33 + //if the LoadHeader function is not zero, the header is invalid
  34 + if (error != 0)
  35 + {
  36 + fid.isValid = 0;
  37 + return fid;
  38 + }
  39 +
  40 + }
  41 +
  42 + //default header to be used in case of writing a file
  43 + else
  44 + {
  45 +
  46 + fid.header.datatype= 4;
  47 + fid.header.nF= nBands;
  48 + fid.header.sX=0;
  49 + fid.header.sY=1;
  50 + fid.header.headersize=0;
  51 + }
  52 + //load the binary file associated with this header
  53 + //get the binary file name
  54 +
  55 + //find the binary file name
  56 +
  57 + //locate the position of the extension
  58 + int pos = filename.find(".hdr");
  59 + string binaryname = filename.substr(0, pos);
  60 + //load the binary file using fopen()
  61 + //store the resulting FILE pointer in the EnviFidClass variable
  62 + mode += "b";
  63 + fid.file = fopen (binaryname.c_str(), mode.c_str());
  64 + if(fid.file == NULL)
  65 + {
  66 + cout<<"Error opening binary file."<<endl;
  67 + fid.isValid = 0;
  68 + return fid;
  69 + }
  70 +
  71 + //set the file name
  72 + fid.fileName = binaryname;
  73 +
  74 + fid.isValid = 1;
  75 + return fid;
  76 +}
  77 +
  78 +static EnviFidClass EnviOpen(EnviFidClass fid, char* mask = NULL, string mode = "r")
  79 +{
  80 + string filename = fid.fileName;
  81 + filename += ".hdr";
  82 +
  83 + return EnviOpen(filename, mask, mode);
  84 +}
  85 +
  86 +#endif
  87 +
... ...
envi/EnviRead.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef ENVIREAD_H
  2 +#define ENVIREAD_H
  3 +
  4 +#include <string>
  5 +#include <iostream>
  6 +#include <fstream>
  7 +#include <ctime>
  8 +#include <sstream>
  9 +#include <stdio.h>
  10 +#include "EnviHeader.h"
  11 +#include "EnviFid.h"
  12 +#include "EnviOpen.h"
  13 +using namespace std;
  14 +
  15 +
  16 +static int EnviRead(float* ptr, int nPixels, EnviFidClass &fid)
  17 +{
  18 + //number of values to read from the file
  19 + int totalPixels = fid.header.sX * fid.header.sY;
  20 +
  21 + //ERROR TEST: make sure that fid.file is a valid file pointer
  22 + if (fid.file == NULL)
  23 + {
  24 + cout<<"Binary file isn't open."<<endl;
  25 + return 0;
  26 + }
  27 +
  28 + //---------MASK implementation
  29 +
  30 +
  31 + int i=0;
  32 + unsigned int sizeRead;
  33 + //iterate through every pixel in the image
  34 + //while we have not encountered the end of the file AND we have not read nPixels pixels
  35 + while(fid.filePos < totalPixels && i < nPixels)
  36 + {
  37 + //if there is no mask OR the current pixel is valid, read it
  38 + if(fid.mask == NULL || fid.mask[fid.filePos] == 1)
  39 + {
  40 + //read the pixel directly from disk
  41 + sizeRead = fread(&ptr[fid.header.nF * i], sizeof(float), fid.header.nF, fid.file);
  42 +
  43 + //catch an error if one occurs
  44 + if (sizeRead != fid.header.nF)
  45 + {
  46 + cout<<"Error during read # "<<i<<endl;
  47 + return 0;
  48 + }
  49 + i++;
  50 + }
  51 +
  52 +
  53 + //otherwise, skip it
  54 + else
  55 + fseek(fid.file, sizeof(float) * fid.header.nF, SEEK_CUR);
  56 +
  57 + //increment the file position (used to index into the mask)
  58 + fid.filePos++;
  59 +
  60 + }
  61 +
  62 + return i;
  63 +}
  64 +
  65 +#endif
... ...
envi/EnviWrite.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef ENVIWRITE_H
  2 +#define ENVIWRITE_H
  3 +#include <fstream>
  4 +#include <stdio.h>
  5 +#include "EnviHeader.h"
  6 +#include "EnviFid.h"
  7 +#include "EnviOpen.h"
  8 +
  9 +// Take a memory pointer, and copy all of that data to disk.
  10 +static int EnviWrite(float* ptr, int npixels, EnviFidClass &fid)
  11 +{
  12 + if(fid.mode == "r") return 0;
  13 +
  14 + //open a binary file for writing
  15 + if (fid.file == NULL)
  16 + {
  17 + cout<<"Binary file isn't open."<<endl;
  18 + return 1;
  19 + }
  20 +
  21 + //allocate an array of zeros equal to the size of a single pixel
  22 + float* zeroPtr = (float*) malloc(sizeof(float) * fid.header.nF);
  23 + memset(zeroPtr, 0, sizeof(float) * fid.header.nF);
  24 +
  25 +
  26 + int j = 0;
  27 + //for each pixel send to the EnviRead function
  28 + for(int p = 0; p<npixels; p++)
  29 + {
  30 +
  31 + //while the mask is zero
  32 + while(fid.mask[fid.filePos] == 0)
  33 + {
  34 + //increment filePos
  35 + fid.filePos++;
  36 + fid.header.sX++;
  37 + j++;
  38 +
  39 + //write a sequence of zeros
  40 + size_t results= fwrite(zeroPtr, sizeof(float), fid.header.nF, fid.file);
  41 + //check for errors
  42 + if (results != fid.header.nF)
  43 + {
  44 + cout<<"Error during read # "<<j<<endl;
  45 + return 0;
  46 + }
  47 + }
  48 +
  49 + //write a pixel
  50 + size_t results= fwrite(&ptr[fid.header.nF * p], sizeof(float), fid.header.nF, fid.file);
  51 + if (results != fid.header.nF)
  52 + {
  53 + cout<<"Error during read # "<<j<<endl;
  54 + return 0;
  55 + }
  56 +
  57 + //increment filePos
  58 + fid.filePos++;
  59 + fid.header.sX++;
  60 + j++;
  61 + }
  62 +
  63 + return j;
  64 +
  65 +}
  66 +
  67 +
  68 +
  69 +
  70 +
  71 +#endif
... ...
envi/rtsEnvi.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef RTS_ENVI_H
  2 +#define RTS_ENVI_H
  3 +
  4 +#include <string>
  5 +#include <iostream>
  6 +#include <fstream>
  7 +#include <sstream>
  8 +
  9 +/* This is a class for reading and writing ENVI binary files. This class will be updated as needed.
  10 +
  11 +What this class CAN currently do:
  12 + *) Write band images to a BSQ file.
  13 + *) Append band images to a BSQ file.
  14 +
  15 +*/
  16 +
  17 +//information from an ENVI header file
  18 +//A good resource can be found here: http://www.exelisvis.com/docs/enviheaderfiles.html
  19 +struct EnviHeader
  20 +{
  21 + std::string name;
  22 +
  23 + std::string description;
  24 +
  25 + unsigned int samples; //x-axis
  26 + unsigned int lines; //y-axis
  27 + unsigned int bands; //spectral axis
  28 + unsigned int header_offset; //header offset for binary file (in bytes)
  29 + std::string file_type; //should be "ENVI Standard"
  30 + unsigned int data_type; //data representation; common value is 4 (32-bit float)
  31 +
  32 + enum interleaveType {BIP, BIL, BSQ}; //bip = Z,X,Y; bil = X,Z,Y; bsq = X,Y,Z
  33 + interleaveType interleave;
  34 +
  35 + std::string sensor_type; //not really used
  36 +
  37 + enum endianType {endianLittle, endianBig};
  38 + endianType byte_order; //little = least significant bit first (most systems)
  39 +
  40 + double x_start, y_start; //coordinates of the upper-left corner of the image
  41 + std::string wavelength_units; //stored wavelength units
  42 + std::string z_plot_titles[2];
  43 +
  44 + double pixel_size[2]; //pixel size along X and Y
  45 +
  46 + std::vector<std::string> band_names; //name for each band in the image
  47 + std::vector<double> wavelength; //wavelength for each band
  48 +
  49 + EnviHeader()
  50 + {
  51 + name = "";
  52 +
  53 + //specify default values for a new or empty ENVI file
  54 + samples = 0;
  55 + lines = 0;
  56 + bands = 0;
  57 + header_offset = 0;
  58 + data_type = 4;
  59 + interleave = BSQ;
  60 + byte_order = endianLittle;
  61 + x_start = y_start = 0;
  62 + pixel_size[0] = pixel_size[1] = 1;
  63 +
  64 + //strings
  65 + file_type = "ENVI Standard";
  66 + sensor_type = "Unknown";
  67 + wavelength_units = "Unknown";
  68 + z_plot_titles[0] = z_plot_titles[1] = "Unknown";
  69 + }
  70 +
  71 + std::string trim(std::string line)
  72 + {
  73 + //trims whitespace from the beginning and end of line
  74 + int start_i, end_i;
  75 + for(start_i=0; start_i < line.length(); start_i++)
  76 + if(line[start_i] != 32)
  77 + {
  78 + break;
  79 + }
  80 +
  81 + for(end_i = line.length()-1; end_i >= start_i; end_i--)
  82 + if(line[end_i] != ' ')
  83 + {
  84 + break;
  85 + }
  86 +
  87 + return line.substr(start_i, end_i - start_i+1);
  88 + }
  89 +
  90 +
  91 + std::string get_token(std::string line)
  92 + {
  93 + //returns a variable name; in this case we look for the '=' sign
  94 + size_t i = line.find_first_of('=');
  95 +
  96 + std::string result;
  97 + if(i != std::string::npos)
  98 + result = trim(line.substr(0, i-1));
  99 +
  100 + return result;
  101 + }
  102 +
  103 + std::string get_data_str(std::string line)
  104 + {
  105 + size_t i = line.find_first_of('=');
  106 +
  107 + std::string result;
  108 + if(i != std::string::npos)
  109 + result = trim(line.substr(i+1));
  110 + else
  111 + {
  112 + std::cout<<"ENVI Header error - data not found for token: "<<get_token(line)<<std::endl;
  113 + exit(1);
  114 + }
  115 + return result;
  116 + }
  117 +
  118 + std::string get_brace_str(std::string token, std::string line, std::ifstream &file)
  119 + {
  120 + //this function assembles all of the characters between curly braces
  121 + //this is how strings are defined within an ENVI file
  122 +
  123 + std::string result;
  124 +
  125 + //first, find the starting brace
  126 + size_t i;
  127 + do
  128 + {
  129 + i = line.find_first_of('{');
  130 + if(i != std::string::npos)
  131 + break;
  132 + }while(file);
  133 +
  134 + //if i is still npos, we have reached the end of the file without a brace...something is wrong
  135 + if(i == std::string::npos)
  136 + {
  137 + std::cout<<"ENVI Header error - string token declared without being defined: "<<token<<std::endl;
  138 + exit(1);
  139 + }
  140 + line = line.substr(i+1);
  141 +
  142 + //copy character data into the result string until we find a closing brace
  143 + while(file)
  144 + {
  145 + i = line.find_first_of('}');
  146 +
  147 +
  148 + if(i != std::string::npos)
  149 + {
  150 + result += line.substr(0, i);
  151 + break;
  152 + }
  153 + else
  154 + result += line;
  155 +
  156 + getline(file, line);
  157 + }
  158 +
  159 + if(i == std::string::npos)
  160 + {
  161 + std::cout<<"ENVI Header error - string token declared without a terminating '}': "<<token<<std::endl;
  162 + exit(1);
  163 + }
  164 +
  165 + return trim(result);
  166 + }
  167 +
  168 + std::vector<std::string> get_string_seq(std::string token, std::string sequence)
  169 + {
  170 + //this function returns a sequence of comma-delimited strings
  171 + std::vector<std::string> result;
  172 +
  173 + std::string entry;
  174 + size_t i;
  175 + do
  176 + {
  177 + i = sequence.find_first_of(',');
  178 + entry = sequence.substr(0, i);
  179 + sequence = sequence.substr(i+1);
  180 + result.push_back(trim(entry));
  181 + }while(i != std::string::npos);
  182 +
  183 + return result;
  184 + }
  185 +
  186 + std::vector<double> get_double_seq(std::string token, std::string sequence)
  187 + {
  188 + //this function returns a sequence of comma-delimited strings
  189 + std::vector<double> result;
  190 +
  191 + std::string entry;
  192 + size_t i;
  193 + do
  194 + {
  195 + i = sequence.find_first_of(',');
  196 + entry = sequence.substr(0, i);
  197 + sequence = sequence.substr(i+1);
  198 + result.push_back(atof(entry.c_str()));
  199 + }while(i != std::string::npos);
  200 +
  201 + return result;
  202 + }
  203 +
  204 + void load(std::string filename, std::string mode = "r")
  205 + {
  206 + name = filename;
  207 +
  208 +
  209 + //open the header file
  210 + std::ifstream file(name.c_str());
  211 +
  212 + if(!file)
  213 + {
  214 + //if we are in read mode, throw an exception
  215 + if(mode == "r")
  216 + {
  217 + std::cout<<"ENVI header file not found: "<<name<<std::endl;
  218 + exit(1);
  219 + }
  220 + else
  221 + {
  222 + //there is no file, so just return the default structure
  223 + return;
  224 + }
  225 + }
  226 +
  227 + //the first line should just be "ENVI"
  228 + std::string line;
  229 + getline(file, line);
  230 + if(line != "ENVI")
  231 + {
  232 + std::cout<<"The header doesn't appear to be an ENVI file. The first line should be 'ENVI'."<<std::endl;
  233 + exit(1);
  234 + }
  235 +
  236 + //for each line in the file, get the token
  237 + std::string token;
  238 + while(file)
  239 + {
  240 +
  241 + //get a line
  242 + getline(file, line);
  243 +
  244 + //get the token
  245 + token = get_token(line);
  246 +
  247 + if(token == "description")
  248 + description = get_brace_str(token, line, file);
  249 + else if(token == "band names")
  250 + {
  251 + std::string string_sequence = get_brace_str(token, line, file);
  252 + band_names = get_string_seq(token, string_sequence);
  253 + }
  254 + else if(token == "wavelength")
  255 + {
  256 + std::string string_sequence = get_brace_str(token, line, file);
  257 + wavelength = get_double_seq(token, string_sequence);
  258 + }
  259 + else if(token == "pixel size")
  260 + {
  261 + std::string string_sequence = get_brace_str(token, line, file);
  262 + std::vector<double> pxsize = get_double_seq(token, string_sequence);
  263 + pixel_size[0] = pxsize[0];
  264 + pixel_size[1] = pxsize[1];
  265 + }
  266 + else if(token == "z plot titles")
  267 + {
  268 + std::string string_sequence = get_brace_str(token, line, file);
  269 + std::vector<std::string> titles = get_string_seq(token, string_sequence);
  270 + z_plot_titles[0] = titles[0];
  271 + z_plot_titles[1] = titles[1];
  272 + }
  273 +
  274 + else if(token == "samples")
  275 + samples = atoi(get_data_str(line).c_str());
  276 + else if(token == "lines")
  277 + {
  278 + lines = atoi(get_data_str(line).c_str());
  279 + std::cout<<"Number of lines loaded: "<<lines<<std::endl;
  280 + }
  281 + else if(token == "bands")
  282 + bands = atoi(get_data_str(line).c_str());
  283 + else if(token == "header offset")
  284 + header_offset = atoi(get_data_str(line).c_str());
  285 + else if(token == "file type")
  286 + file_type = get_data_str(line);
  287 + else if(token == "data type")
  288 + data_type = atoi(get_data_str(line).c_str());
  289 + else if(token == "interleave")
  290 + {
  291 + std::string interleave_str = get_data_str(line);
  292 + if(interleave_str == "bip")
  293 + interleave = BIP;
  294 + else if(interleave_str == "bil")
  295 + interleave = BIL;
  296 + else if(interleave_str == "bsq")
  297 + interleave = BSQ;
  298 + }
  299 + else if(token == "sensor type")
  300 + sensor_type = get_data_str(line);
  301 + else if(token == "byte order")
  302 + byte_order = (endianType)atoi(get_data_str(line).c_str());
  303 + else if(token == "x start")
  304 + x_start = atof(get_data_str(line).c_str());
  305 + else if(token == "y start")
  306 + y_start = atof(get_data_str(line).c_str());
  307 + else if(token == "wavelength units")
  308 + wavelength_units = get_data_str(line);
  309 + }
  310 +
  311 + //close the file
  312 + file.close();
  313 + }
  314 +
  315 + void save(std::string filename)
  316 + {
  317 + //open a file
  318 + std::ofstream outfile(filename.c_str());
  319 +
  320 + //write the ENVI type identifier
  321 + outfile<<"ENVI"<<std::endl;
  322 +
  323 + //output all of the data
  324 + outfile<<"description = {"<<std::endl;
  325 + outfile<<" "<<description<<"}"<<std::endl;
  326 +
  327 + outfile<<"samples = "<<samples<<std::endl;
  328 + outfile<<"lines = "<<lines<<std::endl;
  329 + outfile<<"bands = "<<bands<<std::endl;
  330 + outfile<<"header offset = "<<header_offset<<std::endl;
  331 + outfile<<"file type = "<<file_type<<std::endl;
  332 + outfile<<"data type = "<<data_type<<std::endl;
  333 + outfile<<"interleave = ";
  334 + if(interleave == BIP)
  335 + outfile<<"bip";
  336 + if(interleave == BIL)
  337 + outfile<<"bil";
  338 + if(interleave == BSQ)
  339 + outfile<<"bsq";
  340 + outfile<<std::endl;
  341 + outfile<<"sensor type = "<<sensor_type<<std::endl;
  342 + outfile<<"byte order = "<<byte_order<<std::endl;
  343 + outfile<<"x start = "<<x_start<<std::endl;
  344 + outfile<<"y start = "<<y_start<<std::endl;
  345 + outfile<<"wavelength units = "<<wavelength_units<<std::endl;
  346 + outfile<<"z plot titles = {";
  347 + outfile<<z_plot_titles[0]<<", "<<z_plot_titles[1]<<"}"<<std::endl;
  348 + outfile<<"pixel size = {"<<pixel_size[0]<<", "<<pixel_size[1]<<", units=Meters}"<<std::endl;
  349 + if(band_names.size() > 0)
  350 + {
  351 + outfile<<"band names = {"<<std::endl;
  352 + for(int i=0; i<band_names.size(); i++)
  353 + {
  354 + outfile<<band_names[i];
  355 + if(i < band_names.size() - 1)
  356 + outfile<<", ";
  357 + }
  358 + outfile<<"}"<<std::endl;
  359 + }
  360 + outfile<<"wavelength = {"<<std::endl;
  361 + for(int i=0; i<wavelength.size()-1; i++)
  362 + outfile<<wavelength[i]<<", ";
  363 + outfile<<wavelength.back()<<"}"<<std::endl;
  364 +
  365 + outfile.close();
  366 + }
  367 +
  368 + void save()
  369 + {
  370 + //std::cout<<"ENVI Header Name: "<<name<<std::endl;
  371 + save(name);
  372 + }
  373 +
  374 + //returns the size of a single value (in bytes)
  375 + unsigned int valsize()
  376 + {
  377 + switch(data_type)
  378 + {
  379 + case 1: //1 = 8-bit byte
  380 + return 1;
  381 + case 2: //16-bit signed integer
  382 + case 12: //16-bit unsigned integer
  383 + return 2;
  384 + case 3: //32-bit signed long integer
  385 + case 4: //32-bit floating point
  386 + case 13: //32-bit unsigned long integer
  387 + return 4;
  388 + case 5: //64-bit double-precision floating point
  389 + case 6: //32-bit complex value
  390 + case 14: //64-bit signed long integer
  391 + case 15: //64-bit unsigned long integer
  392 + return 8;
  393 + case 9: //64-bit complex value
  394 + return 16;
  395 + default:
  396 + return 0;
  397 + }
  398 +
  399 + }
  400 +}; //end EnviHeader
  401 +
  402 +class EnviFile
  403 +{
  404 + EnviHeader header;
  405 + FILE* data;
  406 + std::string mode;
  407 +
  408 + //a locked file has alread had data written..certain things can't be changed
  409 + // accessor functions deal with these issues
  410 + bool locked;
  411 +
  412 + void init()
  413 + {
  414 + locked = false;
  415 + data = NULL;
  416 + }
  417 +
  418 + void readonly()
  419 + {
  420 + if(mode == "r")
  421 + {
  422 + std::cout<<"ENVI Error: changes cannot be made to a read-only file."<<std::endl;
  423 + exit(1);
  424 + }
  425 + }
  426 +
  427 + //this function determines if the current software is capable of the requested change
  428 + bool caps(EnviHeader new_header)
  429 + {
  430 + readonly(); //if the file is read-only, throw an exception
  431 +
  432 + //we currently only support BSQ
  433 + if(new_header.interleave == EnviHeader::BIP || new_header.interleave == EnviHeader::BIL)
  434 + {
  435 + std::cout<<"This software only supports BSQ.";
  436 + return false;
  437 + }
  438 +
  439 + //if the number of bands has changed
  440 + if(header.bands != new_header.bands)
  441 + {
  442 + //the new header must be a BSQ file
  443 + if(new_header.interleave != EnviHeader::BSQ)
  444 + {
  445 + std::cout<<"ENVI Error: can only add bands to a BSQ file."<<std::endl;
  446 + return false;
  447 + }
  448 + }
  449 +
  450 + //disallow anything that can't be done when the file is locked
  451 + if(locked)
  452 + {
  453 + if(header.lines != new_header.lines)
  454 + {
  455 + std::cout<<"ENVI Error: cannot change the number of lines from "<<header.lines<<" to "<<new_header.lines<<" in a locked BSQ file."<<std::endl;
  456 + return false;
  457 + }
  458 + if(header.samples != new_header.samples)
  459 + {
  460 + std::cout<<"ENVI Error: cannot change the number of samples in a locked BSQ file."<<std::endl;
  461 + return false;
  462 + }
  463 + if(header.data_type != new_header.data_type)
  464 + {
  465 + std::cout<<"ENVI Error: cannot change the data type of a locked file."<<std::endl;
  466 + return false;
  467 + }
  468 + if(header.byte_order != new_header.byte_order)
  469 + {
  470 + std::cout<<"ENVI Error: cannot change the byte order of a locked file."<<std::endl;
  471 + return false;
  472 + }
  473 + }
  474 +
  475 + return true;
  476 +
  477 + }
  478 +
  479 + EnviHeader load_header(std::string header_file, std::string file_mode)
  480 + {
  481 + EnviHeader new_header;
  482 + new_header.load(header_file, file_mode);
  483 + return new_header;
  484 + }
  485 +
  486 + void set_header(EnviHeader new_header)
  487 + {
  488 + if(caps(new_header))
  489 + header = new_header;
  490 + else
  491 + exit(1);
  492 + }
  493 +
  494 + bool test_exist(std::string filename)
  495 + {
  496 + //tests for the existance of the specified binary file
  497 + FILE* f;
  498 + long sz = 0;
  499 + f = fopen(filename.c_str(), "rb");
  500 + if(f == NULL)
  501 + {
  502 + std::cout<<"ENVI Error: error testing file existance."<<std::endl;
  503 + exit(1);
  504 + }
  505 +
  506 + fseek(f, 0, SEEK_END);
  507 + sz = ftell(f);
  508 + fclose(f);
  509 +
  510 + if(sz>0)
  511 + return true;
  512 + else
  513 + return false;
  514 +
  515 + }
  516 +
  517 +public:
  518 +
  519 + EnviFile()
  520 + {
  521 + init();
  522 + }
  523 +
  524 + EnviFile(std::string filename, std::string file_mode = "r")
  525 + {
  526 + init();
  527 + open(filename, filename + ".hdr", file_mode);
  528 + }
  529 +
  530 + void open(std::string file_name, std::string header_name, std::string file_mode = "r")
  531 + {
  532 + mode = file_mode;
  533 +
  534 + header.name = header_name;
  535 +
  536 + //lock the file if we are loading or appending
  537 + if(mode == "r" || mode == "a" || mode == "r+" || mode == "a+")
  538 + {
  539 + //load the ENVI header - failure will cause an exit
  540 + set_header(load_header(header_name, file_mode));
  541 +
  542 + if(mode == "r" || test_exist(file_name))
  543 + {
  544 + locked = true;
  545 + }
  546 +
  547 + }
  548 + else
  549 + {
  550 + header.name = header_name;
  551 + }
  552 +
  553 + //open the data file
  554 + data = fopen(file_name.c_str(), (mode + "b").c_str());
  555 +
  556 +
  557 + }
  558 +
  559 + void setDescription(std::string desc)
  560 + {
  561 + readonly(); //if the file is read-only, throw an exception
  562 + header.description = desc;
  563 + }
  564 +
  565 + void setZPlotTitles(std::string spectrum, std::string value)
  566 + {
  567 + readonly(); //if the file is read-only, throw an exception
  568 + header.z_plot_titles[0] = spectrum;
  569 + header.z_plot_titles[1] = value;
  570 + }
  571 +
  572 + void setPixelSize(double sx, double sy)
  573 + {
  574 + readonly(); //if the file is read-only, throw an exception
  575 + header.pixel_size[0] = sx;
  576 + header.pixel_size[1] = sy;
  577 + }
  578 +
  579 + void setWavelengthUnits(std::string units)
  580 + {
  581 + readonly(); //if the file is read-only, throw an exception
  582 + header.wavelength_units = units;
  583 + }
  584 +
  585 + void setCoordinates(double x, double y)
  586 + {
  587 + readonly(); //if the file is read-only, throw an exception
  588 + header.x_start = x;
  589 + header.y_start = y;
  590 + }
  591 +
  592 + //FUNCTIONS THAT MAY BE DISALLOWED IN A LOCKED FILE
  593 + void setSamples(unsigned int samples)
  594 + {
  595 + EnviHeader newHeader = header;
  596 + newHeader.samples = samples;
  597 + set_header(newHeader);
  598 + }
  599 +
  600 + void setLines(unsigned int lines)
  601 + {
  602 + EnviHeader newHeader = header;
  603 + newHeader.lines = lines;
  604 + set_header(newHeader);
  605 + }
  606 +
  607 + void setBands(unsigned int bands)
  608 + {
  609 + EnviHeader newHeader = header;
  610 + newHeader.bands = bands;
  611 + set_header(newHeader);
  612 + }
  613 +
  614 + //DATA MANIPULATION
  615 + void addBand(void* band, unsigned int samples, unsigned int lines, double wavelength, std::string band_name ="")
  616 + {
  617 + //add a band to the file
  618 +
  619 + EnviHeader newHeader = header;
  620 + newHeader.bands++;
  621 + newHeader.samples = samples;
  622 + newHeader.lines = lines;
  623 + newHeader.wavelength.push_back(wavelength);
  624 + if(band_name != "")
  625 + newHeader.band_names.push_back(band_name);
  626 + set_header(newHeader);
  627 +
  628 + //copy the data to the file
  629 + fwrite(band, header.valsize(), samples * lines, data);
  630 +
  631 + //since data was written, the file must be locked
  632 + locked = true;
  633 + }
  634 +
  635 + //close file
  636 + void close()
  637 + {
  638 + //std::cout<<"ENVI File Mode: "<<mode<<std::endl;
  639 +
  640 + //close the data stream
  641 + fclose(data);
  642 +
  643 + //close the header
  644 + if(mode != "r")
  645 + header.save();
  646 + }
  647 +
  648 +}; //end EnviFile
  649 +
  650 +
  651 +#endif
... ...
rts/material.h
  Show/Hide comments   View file @a47a23a
... ... @@ -8,7 +8,7 @@
8 8 #include <complex>
9 9 #include <algorithm>
10 10 #include <sstream>
11   -#include "rts/rtcomplex.h"
  11 +#include "rts/rtsComplex.h"
12 12  
13 13 #define PI 3.14159
14 14  
... ... @@ -49,7 +49,7 @@ namespace rts{
49 49 {
50 50 //wavelength (in microns)
51 51 T lambda;
52   - rtcomplex<T> n;
  52 + rtsComplex<T> n;
53 53 };
54 54  
55 55 template <class T>
... ... @@ -196,6 +196,15 @@ namespace rts{
196 196 dispersion.push_back(ri);
197 197 }
198 198  
  199 + void add(T lambda, rtsComplex<T> n)
  200 + {
  201 + refIndex<T> ri;
  202 + ri.lambda = lambda;
  203 + ri.n = n;
  204 +
  205 + dispersion.push_back(ri);
  206 + }
  207 +
199 208 //comparison function for sorting
200 209 static bool compare(refIndex<T> a, refIndex<T> b)
201 210 {
... ... @@ -234,9 +243,9 @@ namespace rts{
234 243  
235 244 #ifdef FFTW_AVAILABLE
236 245 //allocate memory for the FFT
237   - rtcomplex<T>* Chi2 = (rtcomplex<T>*)fftw_malloc(sizeof(rtcomplex<T>) * N * pf);
238   - rtcomplex<T>* Chi2FFT = (rtcomplex<T>*)fftw_malloc(sizeof(rtcomplex<T>) * N * pf);
239   - rtcomplex<T>* Chi1 = (rtcomplex<T>*)fftw_malloc(sizeof(rtcomplex<T>) * N * pf);
  246 + rtsComplex<T>* Chi2 = (rtsComplex<T>*)fftw_malloc(sizeof(rtsComplex<T>) * N * pf);
  247 + rtsComplex<T>* Chi2FFT = (rtsComplex<T>*)fftw_malloc(sizeof(rtsComplex<T>) * N * pf);
  248 + rtsComplex<T>* Chi1 = (rtsComplex<T>*)fftw_malloc(sizeof(rtsComplex<T>) * N * pf);
240 249  
241 250 //create an FFT plan for the forward and inverse transforms
242 251 fftw_plan planForward, planInverse;
... ... @@ -267,8 +276,8 @@ namespace rts{
267 276 }
268 277  
269 278 //use linear interpolation between the start and end points to pad the spectrum
270   - //rtcomplex<T> nMin = dispersion.back();
271   - //rtcomplex<T> nMax = dispersion.front();
  279 + //rtsComplex<T> nMin = dispersion.back();
  280 + //rtsComplex<T> nMax = dispersion.front();
272 281 T a;
273 282 for(int i=N; i<N*pf; i++)
274 283 {
... ... @@ -282,7 +291,7 @@ namespace rts{
282 291 fftw_execute(planForward);
283 292  
284 293 //perform the Hilbert transform in the Fourier domain
285   - rtcomplex<T> j(0, 1);
  294 + rtsComplex<T> j(0, 1);
286 295 for(int i=0; i<N*pf; i++)
287 296 {
288 297 //if w = 0, set the DC component to zero
... ... @@ -375,12 +384,12 @@ namespace rts{
375 384 n0 = n;
376 385 }
377 386  
378   - material(std::string filename, std::string format, T scaleA = 1.0)
  387 + material(std::string filename, std::string format = "", T scaleA = 1.0)
379 388 {
380 389 fromFile(filename, format);
381 390 }
382 391  
383   - void fromFile(std::string filename, std::string format, T scaleA = 1.0)
  392 + void fromFile(std::string filename, std::string format = "", T scaleA = 1.0)
384 393 {
385 394 //clear any previous values
386 395 dispersion.clear();
... ... @@ -402,7 +411,7 @@ namespace rts{
402 411  
403 412 }
404 413  
405   - void fromStr(std::string str, std::string format, T scaleA = 1.0)
  414 + void fromStr(std::string str, std::string format = "", T scaleA = 1.0)
406 415 {
407 416 //create a string stream to process the input data
408 417 std::stringstream ss(str);
... ... @@ -410,9 +419,38 @@ namespace rts{
410 419 //this string will be read line-by-line (where each line is an entry)
411 420 std::string line;
412 421  
413   - //create an entry structure (for now just a basic one)
  422 + //if the format is not provided, see if it is in the file, otherwise use a default
  423 + if(format == "")
  424 + {
  425 + //set a default format of "lambda,n,k"
  426 + format = "microns,n,k";
  427 +
  428 + //see if the first line is a comment
  429 + char c = ss.peek();
  430 + if(c == '#')
  431 + {
  432 + //get the first line
  433 + getline(ss, line);
  434 + //look for a bracket, denoting the start of a format string
  435 + int istart = line.find('[');
  436 + if(istart != std::string::npos)
  437 + {
  438 + //look for a bracket terminating the format string
  439 + int iend = line.find(']');
  440 + if(iend != std::string::npos)
  441 + {
  442 + //read the string between the brackets
  443 + format = line.substr(istart+1, iend - istart - 1);
  444 + }
  445 + }
  446 + }
  447 +
  448 + }
  449 +
414 450 entryType<T> entry(format);
415 451  
  452 + std::cout<<"Loading material with format: "<<format<<std::endl;
  453 +
416 454 T lambda, n, k;
417 455 while(!ss.eof())
418 456 {
... ... @@ -434,24 +472,37 @@ namespace rts{
434 472 }
435 473  
436 474 //convert the material to a string
437   - std::string toStr(std::string format = "microns,n,k")
  475 + std::string toStr(std::string format = "microns,n,k", bool reverse_order = false)
438 476 {
439 477 std::stringstream ss;
440 478 entryType<T> entry(format);
441   - for(unsigned int l=0; l<dispersion.size(); l++)
442   - {
443   - if(l > 0) ss<<std::endl;
444   - ss<<entry.outputEntry(dispersion[l]);
  479 +
  480 + if(reverse_order)
  481 + {
  482 + for(int l=dispersion.size() - 1; l>=0; l--)
  483 + {
  484 + if(l < dispersion.size() - 1) ss<<std::endl;
  485 + ss<<entry.outputEntry(dispersion[l]);
  486 + }
  487 +
  488 + }
  489 + else
  490 + {
  491 + for(unsigned int l=0; l<dispersion.size(); l++)
  492 + {
  493 + if(l > 0) ss<<std::endl;
  494 + ss<<entry.outputEntry(dispersion[l]);
  495 + }
445 496 }
446 497  
447 498 return ss.str();
448 499 }
449 500  
450   - void save(std::string filename, std::string format = "microns,n,k")
  501 + void save(std::string filename, std::string format = "microns,n,k", bool reverse_order = false)
451 502 {
452   - ofstream outfile(filename.c_str());
  503 + std::ofstream outfile(filename.c_str());
453 504 outfile<<"#material file saved as [" + format + "]"<<std::endl;
454   - outfile<<toStr(format)<<std::endl;
  505 + outfile<<toStr(format, reverse_order)<<std::endl;
455 506  
456 507 }
457 508  
... ... @@ -475,7 +526,7 @@ namespace rts{
475 526  
476 527 }
477 528  
478   - rtcomplex<T> getN(T l)
  529 + rtsComplex<T> getN(T l)
479 530 {
480 531 //declare an iterator
481 532 typename std::vector< refIndex<T> >::iterator it;
... ... @@ -499,16 +550,16 @@ namespace rts{
499 550 //std::cout<<lMin<<"----------"<<lMax<<std::endl;
500 551  
501 552 T a = (l - lMin) / (lMax - lMin);
502   - rtcomplex<T> riMin = (*(it - 1)).n;
503   - rtcomplex<T> riMax = (*it).n;
504   - rtcomplex<T> interp;
505   - interp = rtcomplex<T>(a, 0.0) * riMin + rtcomplex<T>(1.0 - a, 0.0) * riMax;
  553 + rtsComplex<T> riMin = (*(it - 1)).n;
  554 + rtsComplex<T> riMax = (*it).n;
  555 + rtsComplex<T> interp;
  556 + interp = rtsComplex<T>(a, 0.0) * riMin + rtsComplex<T>(1.0 - a, 0.0) * riMax;
506 557 return interp;
507 558 }
508 559  
509 560 }
510 561 //interpolate the given lambda value and return the index of refraction
511   - rtcomplex<T> operator()(T l)
  562 + rtsComplex<T> operator()(T l)
512 563 {
513 564 return getN(l);
514 565 }
... ...
rts/rect.h
  Show/Hide comments   View file @a47a23a
... ... @@ -3,8 +3,8 @@
3 3  
4 4 //enable CUDA_CALLABLE macro
5 5 #include "rts/cuda_callable.h"
6   -#include "rts/vector.h"
7   -#include "rts/point.h"
  6 +#include "rts/rtsVector.h"
  7 +#include "rts/rtsPoint.h"
8 8 #include <iostream>
9 9  
10 10 namespace rts{
... ... @@ -27,9 +27,9 @@ struct rect
27 27 T B[N];
28 28 T C[N];*/
29 29  
30   - rts::point<T, N> A;
31   - rts::vector<T, N> X;
32   - rts::vector<T, N> Y;
  30 + rts::rtsPoint<T, N> A;
  31 + rts::rtsVector<T, N> X;
  32 + rts::rtsVector<T, N> Y;
33 33  
34 34  
35 35 CUDA_CALLABLE rect()
... ... @@ -37,7 +37,7 @@ struct rect
37 37  
38 38 }
39 39  
40   - CUDA_CALLABLE rect(point<T, N> a, point<T, N> b, point<T, N> c)
  40 + CUDA_CALLABLE rect(rtsPoint<T, N> a, rtsPoint<T, N> b, rtsPoint<T, N> c)
41 41 {
42 42  
43 43 A = a;
... ... @@ -46,27 +46,27 @@ struct rect
46 46  
47 47 }
48 48  
49   - CUDA_CALLABLE rect(rts::point<T, N> pMin, rts::point<T, N> pMax, rts::vector<T, N> normal)
  49 + CUDA_CALLABLE rect(rts::rtsPoint<T, N> pMin, rts::rtsPoint<T, N> pMax, rts::rtsVector<T, N> normal)
50 50 {
51 51  
52   - //assign the corner point
  52 + //assign the corner rtsPoint
53 53 A = pMin;
54 54  
55 55 //compute the vector from pMin to pMax
56   - rts::vector<T, 3> v0;
  56 + rts::rtsVector<T, 3> v0;
57 57 v0 = pMax - pMin;
58 58  
59 59 //compute the cross product of A and the plane normal
60   - rts::vector<T, 3> v1;
  60 + rts::rtsVector<T, 3> v1;
61 61 v1 = v0.cross(normal);
62 62  
63 63  
64   - //calculate point B
65   - rts::point<T, 3> B;
  64 + //calculate rtsPoint B
  65 + rts::rtsPoint<T, 3> B;
66 66 B = A + v0 * 0.5 + v1 * 0.5;
67 67  
68   - //calculate point C
69   - rts::point<T, 3> C;
  68 + //calculate rtsPoint C
  69 + rts::rtsPoint<T, 3> C;
70 70 C = A + v0 * 0.5 - v1 * 0.5;
71 71  
72 72 //calculate X and Y
... ... @@ -78,7 +78,7 @@ struct rect
78 78  
79 79 }
80 80  
81   - /*CUDA_CALLABLE rect(rts::point<T, N> p, rts::vector<T, N> x, rts::vector<T, N> y, T sx, T sy)
  81 + /*CUDA_CALLABLE rect(rts::rtsPoint<T, N> p, rts::vector<T, N> x, rts::vector<T, N> y, T sx, T sy)
82 82 {
83 83 //This constructor creates a rect given a position, orientation, and size
84 84 // p = center position of the rect
... ... @@ -100,16 +100,16 @@ struct rect
100 100  
101 101 }*/
102 102  
103   - CUDA_CALLABLE rts::point<T, N> p(T a, T b)
  103 + CUDA_CALLABLE rts::rtsPoint<T, N> p(T a, T b)
104 104 {
105   - rts::point<T, N> result;
  105 + rts::rtsPoint<T, N> result;
106 106 //given the two parameters a, b = [0 1], returns the position in world space
107 107 result = A + X * a + Y * b;
108 108  
109 109 return result;
110 110 }
111 111  
112   - CUDA_CALLABLE rts::point<T, N> operator()(T a, T b)
  112 + CUDA_CALLABLE rts::rtsPoint<T, N> operator()(T a, T b)
113 113 {
114 114 return p(a, b);
115 115 }
... ... @@ -126,6 +126,23 @@ struct rect
126 126 return ss.str();
127 127  
128 128 }
  129 +
  130 + CUDA_CALLABLE rect<T, N> operator*(T rhs)
  131 + {
  132 + //scales the plane by a scalar value
  133 +
  134 + //compute the center rtsPoint
  135 + rts::rtsPoint<T, N> c = A + X*0.5 + Y*0.5;
  136 +
  137 + //create the new rectangle
  138 + rect<T, N> result;
  139 + result.X = X * rhs;
  140 + result.Y = Y * rhs;
  141 + result.A = c - result.X*0.5 - result.Y*0.5;
  142 +
  143 + return result;
  144 +
  145 + }
129 146 };
130 147  
131 148 } //end namespace rts
... ...
rts/rtcomplex.h renamed to rts/rtsComplex.h
  Show/Hide comments   View file @a47a23a
... ... @@ -15,54 +15,50 @@ namespace rts
15 15 {
16 16  
17 17 template <class T>
18   -struct rtcomplex
  18 +struct rtsComplex
19 19 {
20 20 T r, i;
21 21  
22 22 //default constructor
23   - CUDA_CALLABLE rtcomplex()
  23 + CUDA_CALLABLE rtsComplex()
24 24 {
25 25 r = 0.0;
26 26 i = 0.0;
27 27 }
28 28  
29 29 //access methods
30   - T real()
  30 + CUDA_CALLABLE T real()
31 31 {
32 32 return r;
33 33 }
34 34  
35   - T real(T r_val)
  35 + CUDA_CALLABLE T real(T r_val)
36 36 {
37 37 r = r_val;
38 38 return r_val;
39 39 }
40 40  
41   - T imag()
  41 + CUDA_CALLABLE T imag()
42 42 {
43 43 return i;
44 44 }
45   - T imag(T i_val)
  45 + CUDA_CALLABLE T imag(T i_val)
46 46 {
47 47 i = i_val;
48 48 return i_val;
49 49 }
50 50  
51   -
52   -
53   -
54   -
55 51 //constructor when given real and imaginary values
56   - CUDA_CALLABLE rtcomplex(T r, T i)
  52 + CUDA_CALLABLE rtsComplex(T r, T i)
57 53 {
58 54 this->r = r;
59 55 this->i = i;
60 56 }
61 57  
62 58 //return the current value multiplied by i
63   - CUDA_CALLABLE rtcomplex<T> imul()
  59 + CUDA_CALLABLE rtsComplex<T> imul()
64 60 {
65   - rtcomplex<T> result;
  61 + rtsComplex<T> result;
66 62 result.r = -i;
67 63 result.i = r;
68 64  
... ... @@ -72,70 +68,70 @@ struct rtcomplex
72 68 //ARITHMETIC OPERATORS--------------------
73 69  
74 70 //binary + operator (returns the result of adding two complex values)
75   - CUDA_CALLABLE rtcomplex<T> operator+ (const rtcomplex<T> rhs)
  71 + CUDA_CALLABLE rtsComplex<T> operator+ (const rtsComplex<T> rhs)
76 72 {
77   - rtcomplex<T> result;
  73 + rtsComplex<T> result;
78 74 result.r = r + rhs.r;
79 75 result.i = i + rhs.i;
80 76 return result;
81 77 }
82 78  
83   - CUDA_CALLABLE rtcomplex<T> operator+ (const T rhs)
  79 + CUDA_CALLABLE rtsComplex<T> operator+ (const T rhs)
84 80 {
85   - rtcomplex<T> result;
  81 + rtsComplex<T> result;
86 82 result.r = r + rhs;
87 83 result.i = i;
88 84 return result;
89 85 }
90 86  
91 87 //binary - operator (returns the result of adding two complex values)
92   - CUDA_CALLABLE rtcomplex<T> operator- (const rtcomplex<T> rhs)
  88 + CUDA_CALLABLE rtsComplex<T> operator- (const rtsComplex<T> rhs)
93 89 {
94   - rtcomplex<T> result;
  90 + rtsComplex<T> result;
95 91 result.r = r - rhs.r;
96 92 result.i = i - rhs.i;
97 93 return result;
98 94 }
99 95  
100 96 //binary - operator (returns the result of adding two complex values)
101   - CUDA_CALLABLE rtcomplex<T> operator- (const T rhs)
  97 + CUDA_CALLABLE rtsComplex<T> operator- (const T rhs)
102 98 {
103   - rtcomplex<T> result;
  99 + rtsComplex<T> result;
104 100 result.r = r - rhs;
105 101 result.i = i;
106 102 return result;
107 103 }
108 104  
109 105 //binary MULTIPLICATION operators (returns the result of multiplying complex values)
110   - CUDA_CALLABLE rtcomplex<T> operator* (const rtcomplex<T> rhs)
  106 + CUDA_CALLABLE rtsComplex<T> operator* (const rtsComplex<T> rhs)
111 107 {
112   - rtcomplex<T> result;
  108 + rtsComplex<T> result;
113 109 result.r = r * rhs.r - i * rhs.i;
114 110 result.i = r * rhs.i + i * rhs.r;
115 111 return result;
116 112 }
117   - CUDA_CALLABLE rtcomplex<T> operator* (const T rhs)
  113 + CUDA_CALLABLE rtsComplex<T> operator* (const T rhs)
118 114 {
119   - return rtcomplex<T>(r * rhs, i * rhs);
  115 + return rtsComplex<T>(r * rhs, i * rhs);
120 116 }
121 117  
122 118 //binary DIVISION operators (returns the result of dividing complex values)
123   - CUDA_CALLABLE rtcomplex<T> operator/ (const rtcomplex<T> rhs)
  119 + CUDA_CALLABLE rtsComplex<T> operator/ (const rtsComplex<T> rhs)
124 120 {
125   - rtcomplex<T> result;
  121 + rtsComplex<T> result;
126 122 T denom = rhs.r * rhs.r + rhs.i * rhs.i;
127 123 result.r = (r * rhs.r + i * rhs.i) / denom;
128 124 result.i = (- r * rhs.i + i * rhs.r) / denom;
129 125  
130 126 return result;
131 127 }
132   - CUDA_CALLABLE rtcomplex<T> operator/ (const T rhs)
  128 + CUDA_CALLABLE rtsComplex<T> operator/ (const T rhs)
133 129 {
134   - return rtcomplex<T>(r / rhs, i / rhs);
  130 + return rtsComplex<T>(r / rhs, i / rhs);
135 131 }
136 132  
137 133 //ASSIGNMENT operators-----------------------------------
138   - CUDA_CALLABLE rtcomplex<T> & operator=(const rtcomplex<T> &rhs)
  134 + CUDA_CALLABLE rtsComplex<T> & operator=(const rtsComplex<T> &rhs)
139 135 {
140 136 //check for self-assignment
141 137 if(this != &rhs)
... ... @@ -145,7 +141,7 @@ struct rtcomplex
145 141 }
146 142 return *this;
147 143 }
148   - CUDA_CALLABLE rtcomplex<T> & operator=(const T &rhs)
  144 + CUDA_CALLABLE rtsComplex<T> & operator=(const T &rhs)
149 145 {
150 146 this->r = rhs;
151 147 this->i = 0;
... ... @@ -154,34 +150,34 @@ struct rtcomplex
154 150 }
155 151  
156 152 //arithmetic assignment operators
157   - CUDA_CALLABLE rtcomplex<T> operator+=(const rtcomplex<T> &rhs)
  153 + CUDA_CALLABLE rtsComplex<T> operator+=(const rtsComplex<T> &rhs)
158 154 {
159 155 *this = *this + rhs;
160 156 return *this;
161 157 }
162   - CUDA_CALLABLE rtcomplex<T> operator+=(const T &rhs)
  158 + CUDA_CALLABLE rtsComplex<T> operator+=(const T &rhs)
163 159 {
164 160 *this = *this + rhs;
165 161 return *this;
166 162 }
167 163  
168   - CUDA_CALLABLE rtcomplex<T> operator*=(const rtcomplex<T> &rhs)
  164 + CUDA_CALLABLE rtsComplex<T> operator*=(const rtsComplex<T> &rhs)
169 165 {
170 166 *this = *this * rhs;
171 167 return *this;
172 168 }
173   - CUDA_CALLABLE rtcomplex<T> operator*=(const T &rhs)
  169 + CUDA_CALLABLE rtsComplex<T> operator*=(const T &rhs)
174 170 {
175 171 *this = *this * rhs;
176 172 return *this;
177 173 }
178 174 //divide and assign
179   - CUDA_CALLABLE rtcomplex<T> operator/=(const rtcomplex<T> &rhs)
  175 + CUDA_CALLABLE rtsComplex<T> operator/=(const rtsComplex<T> &rhs)
180 176 {
181 177 *this = *this / rhs;
182 178 return *this;
183 179 }
184   - CUDA_CALLABLE rtcomplex<T> operator/=(const T &rhs)
  180 + CUDA_CALLABLE rtsComplex<T> operator/=(const T &rhs)
185 181 {
186 182 *this = *this / rhs;
187 183 return *this;
... ... @@ -193,9 +189,9 @@ struct rtcomplex
193 189 return std::sqrt(r * r + i * i);
194 190 }
195 191  
196   - CUDA_CALLABLE rtcomplex<T> log()
  192 + CUDA_CALLABLE rtsComplex<T> log()
197 193 {
198   - rtcomplex<T> result;
  194 + rtsComplex<T> result;
199 195 result.r = std::log(std::sqrt(r * r + i * i));
200 196 result.i = std::atan2(i, r);
201 197  
... ... @@ -203,9 +199,9 @@ struct rtcomplex
203 199 return result;
204 200 }
205 201  
206   - CUDA_CALLABLE rtcomplex<T> exp()
  202 + CUDA_CALLABLE rtsComplex<T> exp()
207 203 {
208   - rtcomplex<T> result;
  204 + rtsComplex<T> result;
209 205  
210 206 T e_r = std::exp(r);
211 207 result.r = e_r * std::cos(i);
... ... @@ -220,18 +216,18 @@ struct rtcomplex
220 216 return pow((double)y);
221 217 }*/
222 218  
223   - CUDA_CALLABLE rtcomplex<T> pow(T y)
  219 + CUDA_CALLABLE rtsComplex<T> pow(T y)
224 220 {
225   - rtcomplex<T> result;
  221 + rtsComplex<T> result;
226 222  
227 223 result = log() * y;
228 224  
229 225 return result.exp();
230 226 }
231 227  
232   - CUDA_CALLABLE rtcomplex<T> sqrt()
  228 + CUDA_CALLABLE rtsComplex<T> sqrt()
233 229 {
234   - rtcomplex<T> result;
  230 + rtsComplex<T> result;
235 231  
236 232 //convert to polar coordinates
237 233 T a = std::sqrt(r*r + i*i);
... ... @@ -257,7 +253,7 @@ struct rtcomplex
257 253 }
258 254  
259 255 //COMPARISON operators
260   - CUDA_CALLABLE bool operator==(rtcomplex<T> rhs)
  256 + CUDA_CALLABLE bool operator==(rtsComplex<T> rhs)
261 257 {
262 258 if(r == rhs.r && i == rhs.i)
263 259 return true;
... ... @@ -301,40 +297,42 @@ struct rtcomplex
301 297  
302 298 };
303 299  
  300 +} //end RTS namespace
  301 +
304 302 //addition
305 303 template<typename T>
306   -CUDA_CALLABLE static rtcomplex<T> operator+(const double a, const rtcomplex<T> b)
  304 +CUDA_CALLABLE static rts::rtsComplex<T> operator+(const double a, const rts::rtsComplex<T> b)
307 305 {
308   - return rtcomplex<T>(a + b.r, b.i);
  306 + return rts::rtsComplex<T>(a + b.r, b.i);
309 307 }
310 308  
311 309 //subtraction with a real value
312 310 template<typename T>
313   -CUDA_CALLABLE static rtcomplex<T> operator-(const double a, const rtcomplex<T> b)
  311 +CUDA_CALLABLE static rts::rtsComplex<T> operator-(const double a, const rts::rtsComplex<T> b)
314 312 {
315   - return rtcomplex<T>(a - b.r, -b.i);
  313 + return rts::rtsComplex<T>(a - b.r, -b.i);
316 314 }
317 315  
318 316 //minus sign
319 317 template<typename T>
320   -CUDA_CALLABLE static rtcomplex<T> operator-(const rtcomplex<T> &rhs)
  318 +CUDA_CALLABLE static rts::rtsComplex<T> operator-(const rts::rtsComplex<T> &rhs)
321 319 {
322   - return rtcomplex<T>(-rhs.r, -rhs.i);
  320 + return rts::rtsComplex<T>(-rhs.r, -rhs.i);
323 321 }
324 322  
325 323 //multiply a T value by a complex value
326 324 template<typename T>
327   -CUDA_CALLABLE static rtcomplex<T> operator*(const double a, const rtcomplex<T> b)
  325 +CUDA_CALLABLE static rts::rtsComplex<T> operator*(const double a, const rts::rtsComplex<T> b)
328 326 {
329   - return rtcomplex<T>(a * b.r, a * b.i);
  327 + return rts::rtsComplex<T>(a * b.r, a * b.i);
330 328 }
331 329  
332 330 //divide a T value by a complex value
333 331 template<typename T>
334   -CUDA_CALLABLE static rtcomplex<T> operator/(const double a, const rtcomplex<T> b)
  332 +CUDA_CALLABLE static rts::rtsComplex<T> operator/(const double a, const rts::rtsComplex<T> b)
335 333 {
336 334 //return complex<T>(a * b.r, a * b.i);
337   - rtcomplex<T> result;
  335 + rts::rtsComplex<T> result;
338 336  
339 337 T denom = b.r * b.r + b.i * b.i;
340 338  
... ... @@ -352,56 +350,62 @@ CUDA_CALLABLE static complex&lt;T&gt; pow(complex&lt;T&gt; x, int y)
352 350 }*/
353 351  
354 352 template<typename T>
355   -CUDA_CALLABLE static rtcomplex<T> pow(rtcomplex<T> x, T y)
  353 +CUDA_CALLABLE static rts::rtsComplex<T> pow(rts::rtsComplex<T> x, T y)
356 354 {
357 355 return x.pow(y);
358 356 }
359 357  
360 358 //log function
361 359 template<typename T>
362   -CUDA_CALLABLE static rtcomplex<T> log(rtcomplex<T> x)
  360 +CUDA_CALLABLE static rts::rtsComplex<T> log(rts::rtsComplex<T> x)
363 361 {
364 362 return x.log();
365 363 }
366 364  
367 365 //exp function
368 366 template<typename T>
369   -CUDA_CALLABLE static rtcomplex<T> exp(rtcomplex<T> x)
  367 +CUDA_CALLABLE static rts::rtsComplex<T> exp(rts::rtsComplex<T> x)
370 368 {
371 369 return x.exp();
372 370 }
373 371  
374 372 //sqrt function
375 373 template<typename T>
376   -CUDA_CALLABLE static rtcomplex<T> sqrt(rtcomplex<T> x)
  374 +CUDA_CALLABLE static rts::rtsComplex<T> sqrt(rts::rtsComplex<T> x)
377 375 {
378 376 return x.sqrt();
379 377 }
380 378  
381 379  
382 380 template <typename T>
383   -CUDA_CALLABLE static T abs(rtcomplex<T> a)
  381 +CUDA_CALLABLE static T abs(rts::rtsComplex<T> a)
384 382 {
385 383 return a.abs();
386 384 }
387 385  
388 386 template <typename T>
389   -CUDA_CALLABLE static T real(rtcomplex<T> a)
  387 +CUDA_CALLABLE static T real(rts::rtsComplex<T> a)
390 388 {
391 389 return a.r;
392 390 }
393 391  
  392 +//template <typename T>
  393 +CUDA_CALLABLE static float real(float a)
  394 +{
  395 + return a;
  396 +}
  397 +
394 398 template <typename T>
395   -CUDA_CALLABLE static T imag(rtcomplex<T> a)
  399 +CUDA_CALLABLE static T imag(rts::rtsComplex<T> a)
396 400 {
397 401 return a.i;
398 402 }
399 403  
400 404 //trigonometric functions
401 405 template<class A>
402   -CUDA_CALLABLE rtcomplex<A> sin(const rtcomplex<A> x)
  406 +CUDA_CALLABLE rts::rtsComplex<A> sin(const rts::rtsComplex<A> x)
403 407 {
404   - rtcomplex<A> result;
  408 + rts::rtsComplex<A> result;
405 409 result.r = std::sin(x.r) * std::cosh(x.i);
406 410 result.i = std::cos(x.r) * std::sinh(x.i);
407 411  
... ... @@ -409,9 +413,9 @@ CUDA_CALLABLE rtcomplex&lt;A&gt; sin(const rtcomplex&lt;A&gt; x)
409 413 }
410 414  
411 415 template<class A>
412   -CUDA_CALLABLE rtcomplex<A> cos(const rtcomplex<A> x)
  416 +CUDA_CALLABLE rts::rtsComplex<A> cos(const rts::rtsComplex<A> x)
413 417 {
414   - rtcomplex<A> result;
  418 + rts::rtsComplex<A> result;
415 419 result.r = std::cos(x.r) * std::cosh(x.i);
416 420 result.i = -(std::sin(x.r) * std::sinh(x.i));
417 421  
... ... @@ -419,11 +423,8 @@ CUDA_CALLABLE rtcomplex&lt;A&gt; cos(const rtcomplex&lt;A&gt; x)
419 423 }
420 424  
421 425  
422   -
423   -} //end RTS namespace
424   -
425 426 template<class A>
426   -std::ostream& operator<<(std::ostream& os, rts::rtcomplex<A> x)
  427 +std::ostream& operator<<(std::ostream& os, rts::rtsComplex<A> x)
427 428 {
428 429 os<<x.toStr();
429 430 return os;
... ...
rts/rtsMatrix.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef RTS_MATRIX_H
  2 +#define RTS_MATRIX_H
  3 +
  4 +//#include "rts/vector.h"
  5 +#include <string.h>
  6 +#include <iostream>
  7 +
  8 +namespace rts
  9 +{
  10 +
  11 +template <class T, int N>
  12 +struct rtsMatrix
  13 +{
  14 + //the matrix will be stored in column-major order (compatible with OpenGL)
  15 + T M[N*N];
  16 +
  17 + rtsMatrix()
  18 + {
  19 + for(int r=0; r<N; r++)
  20 + for(int c=0; c<N; c++)
  21 + if(r == c)
  22 + (*this)(r, c) = 1;
  23 + else
  24 + (*this)(r, c) = 0;
  25 + }
  26 +
  27 + T& operator()(int row, int col)
  28 + {
  29 + return M[col * N + row];
  30 + }
  31 +
  32 + rtsMatrix<T, N> operator=(T rhs)
  33 + {
  34 + int Nsq = N*N;
  35 + for(int i=0; i<Nsq; i++)
  36 + M[i] = rhs;
  37 +
  38 + return *this;
  39 + }
  40 +
  41 + /*rtsMatrix<T, N> operator=(rtsMatrix<T, N> rhs)
  42 + {
  43 + for(int i=0; i<N; i++)
  44 + M[i] = rhs.M[i];
  45 +
  46 + return *this;
  47 + }*/
  48 +
  49 + rtsVector<T, N> operator*(rtsVector<T, N> rhs)
  50 + {
  51 + rtsVector<T, N> result;
  52 +
  53 + for(int r=0; r<N; r++)
  54 + for(int c=0; c<N; c++)
  55 + result[r] += (*this)(r, c) * rhs[c];
  56 +
  57 + return result;
  58 + }
  59 +
  60 + std::string toStr()
  61 + {
  62 + std::stringstream ss;
  63 +
  64 + for(int r = 0; r < N; r++)
  65 + {
  66 + ss<<"| ";
  67 + for(int c=0; c<N; c++)
  68 + {
  69 + ss<<(*this)(r, c)<<" ";
  70 + }
  71 + ss<<"|"<<std::endl;
  72 + }
  73 +
  74 + return ss.str();
  75 + }
  76 +
  77 +
  78 +
  79 +
  80 +};
  81 +
  82 +} //end namespace rts
  83 +
  84 +template <typename T, int N>
  85 +std::ostream& operator<<(std::ostream& os, rts::rtsMatrix<T, N> M)
  86 +{
  87 + os<<M.toStr();
  88 + return os;
  89 +}
  90 +
  91 +#endif
... ...
rts/point.h renamed to rts/rtsPoint.h
  Show/Hide comments   View file @a47a23a
1   -#ifndef RTS_POINT_H
2   -#define RTS_POINT_H
  1 +#ifndef RTS_rtsPoint_H
  2 +#define RTS_rtsPoint_H
3 3  
4 4 //#include "rts/vector.h"
5 5 #include <string.h>
... ... @@ -8,17 +8,17 @@ namespace rts
8 8 {
9 9  
10 10 template <class T, int N>
11   -struct point
  11 +struct rtsPoint
12 12 {
13 13 T p[N];
14 14  
15   - CUDA_CALLABLE point()
  15 + CUDA_CALLABLE rtsPoint()
16 16 {
17 17  
18 18 }
19 19  
20 20 //efficiency constructor, makes construction easier for 1D-4D vectors
21   - CUDA_CALLABLE point(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
  21 + CUDA_CALLABLE rtsPoint(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
22 22 {
23 23 if(N >= 1)
24 24 p[0] = x;
... ... @@ -31,48 +31,48 @@ struct point
31 31 }
32 32  
33 33 //arithmetic operators
34   - CUDA_CALLABLE rts::point<T, N> operator+(rts::vector<T, N> v)
  34 + CUDA_CALLABLE rts::rtsPoint<T, N> operator+(rts::rtsVector<T, N> v)
35 35 {
36   - rts::point<T, N> r;
  36 + rts::rtsPoint<T, N> r;
37 37  
38   - //calculate the position of the resulting point
  38 + //calculate the position of the resulting rtsPoint
39 39 for(int i=0; i<N; i++)
40 40 r.p[i] = p[i] + v.v[i];
41 41  
42 42 return r;
43 43 }
44   - CUDA_CALLABLE rts::point<T, N> operator-(rts::vector<T, N> v)
  44 + CUDA_CALLABLE rts::rtsPoint<T, N> operator-(rts::rtsVector<T, N> v)
45 45 {
46   - rts::point<T, N> r;
  46 + rts::rtsPoint<T, N> r;
47 47  
48   - //calculate the position of the resulting point
  48 + //calculate the position of the resulting rtsPoint
49 49 for(int i=0; i<N; i++)
50 50 r.p[i] = p[i] - v.v[i];
51 51  
52 52 return r;
53 53 }
54   - CUDA_CALLABLE rts::vector<T, N> operator-(rts::point<T, N> rhs)
  54 + CUDA_CALLABLE rts::rtsVector<T, N> operator-(rts::rtsPoint<T, N> rhs)
55 55 {
56   - rts::vector<T, N> r;
  56 + rts::rtsVector<T, N> r;
57 57  
58   - //calculate the position of the resulting point
  58 + //calculate the position of the resulting rtsPoint
59 59 for(int i=0; i<N; i++)
60 60 r.v[i] = p[i] - rhs.p[i];
61 61  
62 62 return r;
63 63 }
64   - CUDA_CALLABLE rts::point<T, N> operator*(T rhs)
  64 + CUDA_CALLABLE rts::rtsPoint<T, N> operator*(T rhs)
65 65 {
66   - rts::point<T, N> r;
  66 + rts::rtsPoint<T, N> r;
67 67  
68   - //calculate the position of the resulting point
  68 + //calculate the position of the resulting rtsPoint
69 69 for(int i=0; i<N; i++)
70 70 r.p[i] = p[i] * rhs;
71 71  
72 72 return r;
73 73 }
74 74  
75   - CUDA_CALLABLE point(const T(&data)[N])
  75 + CUDA_CALLABLE rtsPoint(const T(&data)[N])
76 76 {
77 77 memcpy(p, data, sizeof(T) * N);
78 78 }
... ... @@ -104,7 +104,7 @@ struct point
104 104 } //end namespace rts
105 105  
106 106 template <typename T, int N>
107   -std::ostream& operator<<(std::ostream& os, rts::point<T, N> p)
  107 +std::ostream& operator<<(std::ostream& os, rts::rtsPoint<T, N> p)
108 108 {
109 109 os<<p.toStr();
110 110 return os;
... ... @@ -112,9 +112,9 @@ std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, rts::point&lt;T, N&gt; p)
112 112  
113 113 //arithmetic
114 114 template <typename T, int N>
115   -CUDA_CALLABLE rts::point<T, N> operator*(T lhs, rts::point<T, N> rhs)
  115 +CUDA_CALLABLE rts::rtsPoint<T, N> operator*(T lhs, rts::rtsPoint<T, N> rhs)
116 116 {
117   - rts::point<T, N> r;
  117 + rts::rtsPoint<T, N> r;
118 118  
119 119 return rhs * lhs;
120 120 }
... ...
rts/rtsProgressBar.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef RTS_PROGRESSBAR_H
  2 +#define RTS_PROGRESSBAR_H
  3 +
  4 +#include <iostream>
  5 +#include <sstream>
  6 +using namespace std;
  7 +
  8 +static void rtsProgressBar(int percent)
  9 +{
  10 + stringstream bar;
  11 + static int x = 0;
  12 + string slash[4];
  13 + slash[0] = "\\";
  14 + slash[1] = "-";
  15 + slash[2] = "/";
  16 + slash[3] = "|";
  17 + bar<<"[";
  18 + for(int i=0; i<40; i++)
  19 + {
  20 + if(percent > (float)i/(float)40 *100)
  21 + bar << "*";
  22 + else
  23 + bar << " ";
  24 + }
  25 + bar<<"]";
  26 + cout << "\r"; // carriage return back to beginning of line
  27 + cout << bar.str() << " " << slash[x] << " " << percent << " %"; // print the bars and percentage
  28 + x++; // increment to make the slash appear to rotate
  29 + if(x == 4)
  30 + x = 0; // reset slash animation
  31 +}
  32 +
  33 +#endif
... ...
rts/rtsQuad.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#ifndef RTS_RECT_H
  2 +#define RTS_RECT_H
  3 +
  4 +//enable CUDA_CALLABLE macro
  5 +#include "rts/cuda_callable.h"
  6 +#include "rts/rtsVector.h"
  7 +#include "rts/rtsPoint.h"
  8 +#include <iostream>
  9 +
  10 +namespace rts{
  11 +
  12 +//template for a rtsQuadangle class in ND space
  13 +template <class T, int N>
  14 +struct rtsQuad
  15 +{
  16 + /*
  17 + C------------------>O
  18 + ^ ^
  19 + | |
  20 + Y |
  21 + | |
  22 + | |
  23 + A---------X-------->B
  24 + */
  25 +
  26 + /*T A[N];
  27 + T B[N];
  28 + T C[N];*/
  29 +
  30 + rts::rtsPoint<T, N> A;
  31 + rts::rtsVector<T, N> X;
  32 + rts::rtsVector<T, N> Y;
  33 +
  34 +
  35 + CUDA_CALLABLE rtsQuad()
  36 + {
  37 +
  38 + }
  39 +
  40 + CUDA_CALLABLE rtsQuad(rtsPoint<T, N> a, rtsPoint<T, N> b, rtsPoint<T, N> c)
  41 + {
  42 +
  43 + A = a;
  44 + X = b - a;
  45 + Y = c - a;
  46 +
  47 + }
  48 +
  49 + CUDA_CALLABLE rtsQuad(rts::rtsPoint<T, N> pMin, rts::rtsPoint<T, N> pMax, rts::rtsVector<T, N> normal)
  50 + {
  51 +
  52 + //assign the corner rtsPoint
  53 + A = pMin;
  54 +
  55 + //compute the vector from pMin to pMax
  56 + rts::rtsVector<T, 3> v0;
  57 + v0 = pMax - pMin;
  58 +
  59 + //compute the cross product of A and the plane normal
  60 + rts::rtsVector<T, 3> v1;
  61 + v1 = v0.cross(normal);
  62 +
  63 +
  64 + //calculate rtsPoint B
  65 + rts::rtsPoint<T, 3> B;
  66 + B = A + v0 * 0.5 + v1 * 0.5;
  67 +
  68 + //calculate rtsPoint C
  69 + rts::rtsPoint<T, 3> C;
  70 + C = A + v0 * 0.5 - v1 * 0.5;
  71 +
  72 + //calculate X and Y
  73 + X = B - A;
  74 + Y = C - A;
  75 +
  76 +
  77 +
  78 +
  79 + }
  80 +
  81 + /*CUDA_CALLABLE rtsQuad(rts::rtsPoint<T, N> p, rts::vector<T, N> x, rts::vector<T, N> y, T sx, T sy)
  82 + {
  83 + //This constructor creates a rtsQuad given a position, orientation, and size
  84 + // p = center position of the rtsQuad
  85 + // x = x-axis for the rtsQuadangle
  86 + // y = y-axis for the rtsQuadangle
  87 + // sx = size of the rtsQuad along the A-B axis
  88 + // sy = size of the rtsQuad along the A-C axis
  89 +
  90 + //normalize x and y
  91 + rts::vector<T, N> nx = x.norm();
  92 + rts::vector<T, N> ny = y.norm();
  93 +
  94 + //compute X and Y
  95 + X = sx * x;
  96 + Y = sy * y;
  97 +
  98 + //compute A
  99 + A = p - 0.5 * X - 0.5 * Y;
  100 +
  101 + }*/
  102 +
  103 + CUDA_CALLABLE rts::rtsPoint<T, N> p(T a, T b)
  104 + {
  105 + rts::rtsPoint<T, N> result;
  106 + //given the two parameters a, b = [0 1], returns the position in world space
  107 + result = A + X * a + Y * b;
  108 +
  109 + return result;
  110 + }
  111 +
  112 + CUDA_CALLABLE rts::rtsPoint<T, N> operator()(T a, T b)
  113 + {
  114 + return p(a, b);
  115 + }
  116 +
  117 + std::string toStr()
  118 + {
  119 + std::stringstream ss;
  120 +
  121 + ss<<"A = "<<A<<std::endl;
  122 + ss<<"B = "<<A + X<<std::endl;
  123 + ss<<"C = "<<A + X + Y<<std::endl;
  124 + ss<<"D = "<<A + Y<<std::endl;
  125 +
  126 + return ss.str();
  127 +
  128 + }
  129 +
  130 + CUDA_CALLABLE rtsQuad<T, N> operator*(T rhs)
  131 + {
  132 + //scales the plane by a scalar value
  133 +
  134 + //compute the center rtsPoint
  135 + rts::rtsPoint<T, N> c = A + X*0.5 + Y*0.5;
  136 +
  137 + //create the new rtsQuadangle
  138 + rtsQuad<T, N> result;
  139 + result.X = X * rhs;
  140 + result.Y = Y * rhs;
  141 + result.A = c - result.X*0.5 - result.Y*0.5;
  142 +
  143 + return result;
  144 +
  145 + }
  146 +};
  147 +
  148 +} //end namespace rts
  149 +
  150 +template <typename T, int N>
  151 +std::ostream& operator<<(std::ostream& os, rts::rtsQuad<T, N> R)
  152 +{
  153 + os<<R.toStr();
  154 + return os;
  155 +}
  156 +
  157 +
  158 +#endif
0 159 \ No newline at end of file
... ...
rts/rtsQuaternion.h 0 → 100644
  Show/Hide comments   View file @a47a23a
  1 +#include "rtsMatrix.h"
  2 +
  3 +#ifndef RTS_QUATERNION_H
  4 +#define RTS_QUATERNION_H
  5 +
  6 +namespace rts{
  7 +
  8 +template<typename T>
  9 +class rtsQuaternion
  10 +{
  11 +public:
  12 + T w;
  13 + T x;
  14 + T y;
  15 + T z;
  16 +
  17 + void normalize();
  18 + void CreateRotation(T theta, T axis_x, T axis_y, T axis_z);
  19 + void CreateRotation(T theta, rtsVector<T, 3> axis);
  20 + rtsQuaternion<T> operator*(rtsQuaternion<T> &rhs);
  21 + rtsMatrix<T, 3> toMatrix();
  22 +
  23 +
  24 + rtsQuaternion();
  25 + rtsQuaternion(T w, T x, T y, T z);
  26 +
  27 +};
  28 +
  29 +template<typename T>
  30 +void rtsQuaternion<T>::normalize()
  31 +{
  32 + double length=sqrt(w*w + x*x + y*y + z*z);
  33 + w=w/length;
  34 + x=x/length;
  35 + y=y/length;
  36 + z=z/length;
  37 +}
  38 +
  39 +template<typename T>
  40 +void rtsQuaternion<T>::CreateRotation(T theta, T axis_x, T axis_y, T axis_z)
  41 +{
  42 + //assign the given Euler rotation to this quaternion
  43 + w = cos(theta/2.0);
  44 + x = axis_x*sin(theta/2.0);
  45 + y = axis_y*sin(theta/2.0);
  46 + z = axis_z*sin(theta/2.0);
  47 +}
  48 +
  49 +template<typename T>
  50 +void rtsQuaternion<T>::CreateRotation(T theta, rtsVector<T, 3> axis)
  51 +{
  52 + CreateRotation(theta, axis[0], axis[1], axis[2]);
  53 +}
  54 +
  55 +template<typename T>
  56 +rtsQuaternion<T> rtsQuaternion<T>::operator *(rtsQuaternion<T> &param)
  57 +{
  58 + float A, B, C, D, E, F, G, H;
  59 +
  60 +
  61 + A = (w + x)*(param.w + param.x);
  62 + B = (z - y)*(param.y - param.z);
  63 + C = (w - x)*(param.y + param.z);
  64 + D = (y + z)*(param.w - param.x);
  65 + E = (x + z)*(param.x + param.y);
  66 + F = (x - z)*(param.x - param.y);
  67 + G = (w + y)*(param.w - param.z);
  68 + H = (w - y)*(param.w + param.z);
  69 +
  70 + rtsQuaternion<T> result;
  71 + result.w = B + (-E - F + G + H) /2;
  72 + result.x = A - (E + F + G + H)/2;
  73 + result.y = C + (E - F + G - H)/2;
  74 + result.z = D + (E - F - G + H)/2;
  75 +
  76 + return result;
  77 +}
  78 +
  79 +template<typename T>
  80 +rtsMatrix<T, 3> rtsQuaternion<T>::toMatrix()
  81 +{
  82 + rtsMatrix<T, 3> result;
  83 +
  84 +
  85 + double wx, wy, wz, xx, yy, yz, xy, xz, zz, x2, y2, z2;
  86 +
  87 +
  88 + // calculate coefficients
  89 + x2 = x + x; y2 = y + y;
  90 + z2 = z + z;
  91 + xx = x * x2; xy = x * y2; xz = x * z2;
  92 + yy = y * y2; yz = y * z2; zz = z * z2;
  93 + wx = w * x2; wy = w * y2; wz = w * z2;
  94 +
  95 + result(0, 0) = 1.0 - (yy + zz);
  96 + result(0, 1) = xy - wz;
  97 + //m[0][0] = 1.0 - (yy + zz); m[1][0] = xy - wz;
  98 +
  99 + result(0, 2) = xz + wy;
  100 + //result(0, 3) = 0.0;
  101 + //m[2][0] = xz + wy; m[3][0] = 0.0;
  102 +
  103 + result(1, 0) = xy + wz;
  104 + result(1, 1) = 1.0 - (xx + zz);
  105 + //m[0][1] = xy + wz; m[1][1] = 1.0 - (xx + zz);
  106 +
  107 + result(1, 2) = yz - wx;
  108 + //result(1, 3) = 0.0;
  109 + //m[2][1] = yz - wx; m[3][1] = 0.0;
  110 +
  111 + result(2, 0) = xz - wy;
  112 + result(2, 1) = yz + wx;
  113 + //m[0][2] = xz - wy; m[1][2] = yz + wx;
  114 +
  115 + result(2, 2) = 1.0 - (xx + yy);
  116 + //result(3, 2) = 0.0;
  117 + //m[2][2] = 1.0 - (xx + yy); m[3][2] = 0.0;
  118 +
  119 + /*result(3, 0) = 0.0;
  120 + result(3, 1) = 0.0;
  121 + result(3, 2) = 0.0;
  122 + result(3, 3) = 1.0;*/
  123 + //m[0][3] = 0; m[1][3] = 0;
  124 + //m[2][3] = 0; m[3][3] = 1;
  125 + /*
  126 + double* orientationmatrix=(double*)m;
  127 + char c;
  128 +
  129 +
  130 + double* result=new double[16];
  131 + double* array=(double*)m;
  132 + for(int i=0; i<16; i++)
  133 + result[i]=array[i];
  134 + */
  135 +
  136 + return result;
  137 +}
  138 +
  139 +template<typename T>
  140 +rtsQuaternion<T>::rtsQuaternion()
  141 +{
  142 + w=0.0; x=0.0; y=0.0; z=0.0;
  143 +}
  144 +
  145 +template<typename T>
  146 +rtsQuaternion<T>::rtsQuaternion(T c, T i, T j, T k)
  147 +{
  148 + w=c; x=i; y=j; z=k;
  149 +}
  150 +
  151 +} //end rts namespace
  152 +
  153 +#endif
... ...
rts/vector.h renamed to rts/rtsVector.h
  Show/Hide comments   View file @a47a23a
... ... @@ -10,17 +10,19 @@ namespace rts
10 10  
11 11  
12 12 template <class T, int N>
13   -struct vector
  13 +struct rtsVector
14 14 {
15 15 T v[N];
16 16  
17   - CUDA_CALLABLE vector()
  17 + CUDA_CALLABLE rtsVector()
18 18 {
19 19 //memset(v, 0, sizeof(T) * N);
  20 + for(int i=0; i<N; i++)
  21 + v[i] = 0;
20 22 }
21 23  
22 24 //efficiency constructor, makes construction easier for 1D-4D vectors
23   - CUDA_CALLABLE vector(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
  25 + CUDA_CALLABLE rtsVector(T x, T y = (T)0.0, T z = (T)0.0, T w = (T)0.0)
24 26 {
25 27 if(N >= 1)
26 28 v[0] = x;
... ... @@ -32,7 +34,7 @@ struct vector
32 34 v[3] = w;
33 35 }
34 36  
35   - CUDA_CALLABLE vector(const T(&data)[N])
  37 + CUDA_CALLABLE rtsVector(const T(&data)[N])
36 38 {
37 39 memcpy(v, data, sizeof(T) * N);
38 40 }
... ... @@ -49,25 +51,25 @@ struct vector
49 51  
50 52 }
51 53  
52   - CUDA_CALLABLE vector<T, N> cart2sph()
  54 + CUDA_CALLABLE rtsVector<T, N> cart2sph()
53 55 {
54 56 //convert the vector from cartesian to spherical coordinates
55   - //x, y, z -> r, theta, phi
  57 + //x, y, z -> r, theta, phi (where theta = 0 to 2*pi)
56 58  
57   - vector<T, N> sph;
58   - sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
59   - sph[1] = std::atan2(v[1], v[0]);
  59 + rtsVector<T, N> sph;
  60 + sph[0] = std::sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);
  61 + sph[1] = std::atan2(v[1], v[0]);
60 62 sph[2] = std::acos(v[2] / sph[0]);
61 63  
62 64 return sph;
63 65 }
64 66  
65   - CUDA_CALLABLE vector<T, N> sph2cart()
  67 + CUDA_CALLABLE rtsVector<T, N> sph2cart()
66 68 {
67 69 //convert the vector from cartesian to spherical coordinates
68   - //r, theta, phi -> x, y, z
  70 + //r, theta, phi -> x, y, z (where theta = 0 to 2*pi)
69 71  
70   - vector<T, N> cart;
  72 + rtsVector<T, N> cart;
71 73 cart[0] = v[0] * std::cos(v[1]) * std::sin(v[2]);
72 74 cart[1] = v[0] * std::sin(v[1]) * std::sin(v[2]);
73 75 cart[2] = v[0] * std::cos(v[2]);
... ... @@ -75,10 +77,10 @@ struct vector
75 77 return cart;
76 78 }
77 79  
78   - CUDA_CALLABLE vector<T, N> norm()
  80 + CUDA_CALLABLE rtsVector<T, N> norm()
79 81 {
80 82 //compute and return the vector norm
81   - vector<T, N> result;
  83 + rtsVector<T, N> result;
82 84  
83 85 //compute the vector length
84 86 T l = len();
... ... @@ -92,9 +94,9 @@ struct vector
92 94 return result;
93 95 }
94 96  
95   - CUDA_CALLABLE vector<T, 3> cross(vector<T, 3> rhs)
  97 + CUDA_CALLABLE rtsVector<T, 3> cross(rtsVector<T, 3> rhs)
96 98 {
97   - vector<T, 3> result;
  99 + rtsVector<T, 3> result;
98 100  
99 101 //compute the cross product (only valid for 3D vectors)
100 102 result[0] = v[1] * rhs[2] - v[2] * rhs[1];
... ... @@ -103,8 +105,8 @@ struct vector
103 105  
104 106 return result;
105 107 }
106   -
107   - CUDA_CALLABLE T dot(vector<T, N> rhs)
  108 +
  109 + CUDA_CALLABLE T dot(rtsVector<T, N> rhs)
108 110 {
109 111 T result = (T)0;
110 112  
... ... @@ -112,40 +114,40 @@ struct vector
112 114 result += v[i] * rhs.v[i];
113 115  
114 116 return result;
115   -
  117 +
116 118 }
117 119  
118 120 //arithmetic
119   - CUDA_CALLABLE vector<T, N> operator+(vector<T, N> rhs)
  121 + CUDA_CALLABLE rtsVector<T, N> operator+(rtsVector<T, N> rhs)
120 122 {
121   - vector<T, N> result;
  123 + rtsVector<T, N> result;
122 124  
123 125 for(int i=0; i<N; i++)
124 126 result.v[i] = v[i] + rhs.v[i];
125 127  
126 128 return result;
127 129 }
128   - CUDA_CALLABLE vector<T, N> operator-(vector<T, N> rhs)
  130 + CUDA_CALLABLE rtsVector<T, N> operator-(rtsVector<T, N> rhs)
129 131 {
130   - vector<T, N> result;
  132 + rtsVector<T, N> result;
131 133  
132 134 for(int i=0; i<N; i++)
133 135 result.v[i] = v[i] - rhs.v[i];
134 136  
135 137 return result;
136 138 }
137   - CUDA_CALLABLE vector<T, N> operator*(T rhs)
  139 + CUDA_CALLABLE rtsVector<T, N> operator*(T rhs)
138 140 {
139   - vector<T, N> result;
  141 + rtsVector<T, N> result;
140 142  
141 143 for(int i=0; i<N; i++)
142 144 result.v[i] = v[i] * rhs;
143 145  
144 146 return result;
145 147 }
146   - CUDA_CALLABLE vector<T, N> operator/(T rhs)
  148 + CUDA_CALLABLE rtsVector<T, N> operator/(T rhs)
147 149 {
148   - vector<T, N> result;
  150 + rtsVector<T, N> result;
149 151  
150 152 for(int i=0; i<N; i++)
151 153 result.v[i] = v[i] / rhs;
... ... @@ -181,7 +183,7 @@ struct vector
181 183 } //end namespace rts
182 184  
183 185 template <typename T, int N>
184   -std::ostream& operator<<(std::ostream& os, rts::vector<T, N> v)
  186 +std::ostream& operator<<(std::ostream& os, rts::rtsVector<T, N> v)
185 187 {
186 188 os<<v.toStr();
187 189 return os;
... ... @@ -189,9 +191,9 @@ std::ostream&amp; operator&lt;&lt;(std::ostream&amp; os, rts::vector&lt;T, N&gt; v)
189 191  
190 192 //arithmetic operators
191 193 template <typename T, int N>
192   -CUDA_CALLABLE rts::vector<T, N> operator-(rts::vector<T, N> v)
  194 +CUDA_CALLABLE rts::rtsVector<T, N> operator-(rts::rtsVector<T, N> v)
193 195 {
194   - rts::vector<T, N> r;
  196 + rts::rtsVector<T, N> r;
195 197  
196 198 //negate the vector
197 199 for(int i=0; i<N; i++)
... ... @@ -201,9 +203,9 @@ CUDA_CALLABLE rts::vector&lt;T, N&gt; operator-(rts::vector&lt;T, N&gt; v)
201 203 }
202 204  
203 205 template <typename T, int N>
204   -CUDA_CALLABLE rts::vector<T, N> operator*(T lhs, rts::vector<T, N> rhs)
  206 +CUDA_CALLABLE rts::rtsVector<T, N> operator*(T lhs, rts::rtsVector<T, N> rhs)
205 207 {
206   - rts::vector<T, N> r;
  208 + rts::rtsVector<T, N> r;
207 209  
208 210 return rhs * lhs;
209 211 }
... ...
rts/sbessel.h
  Show/Hide comments   View file @a47a23a
... ... @@ -5,10 +5,11 @@
5 5  
6 6 namespace rts{
7 7  
8   -#define RTS_BESSEL_CONVERGENCE_FLOAT 0.00045
  8 +#define RTS_BESSEL_CONVERGENCE_FLOAT 0.0145
  9 +#define RTS_BESSEL_MAXIMUM_FLOAT -1e33
9 10  
10 11 template <typename T>
11   -CUDA_CALLABLE void sbesselj(int n, rtcomplex<T> x, rtcomplex<T>* j)
  12 +CUDA_CALLABLE void sbesselj(int n, rtsComplex<T> x, rtsComplex<T>* j)
12 13 {
13 14 //compute the first bessel function
14 15 if(n >= 0)
... ... @@ -34,7 +35,7 @@ CUDA_CALLABLE void sbesselj(int n, rtcomplex&lt;T&gt; x, rtcomplex&lt;T&gt;* j)
34 35 }
35 36  
36 37 template <typename T>
37   -CUDA_CALLABLE void sbessely(int n, rtcomplex<T> x, rtcomplex<T>* y)
  38 +CUDA_CALLABLE void sbessely(int n, rtsComplex<T> x, rtsComplex<T>* y)
38 39 {
39 40 //compute the first bessel function
40 41 if(n >= 0)
... ... @@ -54,14 +55,14 @@ CUDA_CALLABLE void sbessely(int n, rtcomplex&lt;T&gt; x, rtcomplex&lt;T&gt;* y)
54 55  
55 56 //spherical Hankel functions of the first kind
56 57 template <typename T>
57   -CUDA_CALLABLE void sbesselh1(int n, rtcomplex<T> x, rtcomplex<T>* h)
  58 +CUDA_CALLABLE void sbesselh1(int n, rtsComplex<T> x, rtsComplex<T>* h)
58 59 {
59 60 //compute j_0 and j_1
60   - rtcomplex<T> j[2];
  61 + rtsComplex<T> j[2];
61 62 sbesselj(1, x, j);
62 63  
63 64 //compute y_0 and y_1
64   - rtcomplex<T> y[2];
  65 + rtsComplex<T> y[2];
65 66 sbessely(1, x, y);
66 67  
67 68 //compute the first-order Hhankel function
... ... @@ -82,55 +83,60 @@ CUDA_CALLABLE void sbesselh1(int n, rtcomplex&lt;T&gt; x, rtcomplex&lt;T&gt;* h)
82 83 template <typename T>
83 84 CUDA_CALLABLE void init_sbesselj(T x, T* j)
84 85 {
85   - /*if(x < EPSILON_FLOAT)
86   - {
87   - j[0] = (T)1;
88   - j[1] = (T)0;
89   - return;
90   - }*/
91   -
92 86 //compute the first 2 bessel functions
93   - j[0] = std::sin(x) / x;
94   -
95   - j[1] = j[0] / x - std::cos(x) / x;
  87 + j[0] = sin(x) / x;
96 88  
97   - //deal with the degenerate case of x = 0
98   - /*if(isnan(j[0]))
99   - {
100   - j[0] = (T)1;
101   - j[1] = (T)0;
102   - }*/
103   -
  89 + j[1] = j[0] / x - cos(x) / x;
104 90 }
105 91  
106 92 template <typename T>
107   -CUDA_CALLABLE void shift_sbesselj(int n, T x, T* j)//, T stability = 1.4)
  93 +CUDA_CALLABLE void init_sbessely(T x, T* y)
108 94 {
  95 + //compute the first 2 bessel functions
  96 + y[0] = -cos(x) / x;
109 97  
110   - /*if(x < EPSILON_FLOAT)
111   - {
112   - j[0] = j[1];
113   - j[1] = (T)0;
114   - return;
115   - }*/
  98 + y[1] = y[0] / x - sin(x) / x;
  99 +}
116 100  
117   - T jnew;
  101 +template <typename T>
  102 +CUDA_CALLABLE void shift_sbesselj(int n, T x, T* b)//, T stability = 1.4)
  103 +{
  104 +
  105 + T bnew;
118 106  
119 107 //compute the next (order n) Bessel function
120   - jnew = ((2 * n - 1) * j[1])/x - j[0];
  108 + bnew = ((2 * n - 1) * b[1])/x - b[0];
121 109  
122 110 //if(n > stability*x)
123   - if(n > x)
124   - if(jnew < RTS_BESSEL_CONVERGENCE_FLOAT)
125   - jnew = (T)0;
  111 + if(n > real(x))
  112 + if(real(bnew) < RTS_BESSEL_CONVERGENCE_FLOAT)
  113 + bnew = 0.0;
  114 +
  115 + //shift and add the new value to the array
  116 + b[0] = b[1];
  117 + b[1] = bnew;
  118 +}
  119 +
  120 +template <typename T>
  121 +CUDA_CALLABLE void shift_sbessely(int n, T x, T* b)//, T stability = 1.4)
  122 +{
  123 +
  124 + T bnew;
  125 +
  126 + //compute the next (order n) Bessel function
  127 + bnew = ((2 * n - 1) * b[1])/x - b[0];
  128 +
  129 + if(bnew < RTS_BESSEL_MAXIMUM_FLOAT ||
  130 + (n > x && bnew > 0))
  131 + {
  132 + bnew = (T)0;
  133 + b[1] = (T)0;
  134 + }
126 135  
127   - //deal with the degenerate case of x = 0
128   - //if(isnan(jnew))
129   - // jnew = (T)0;
130 136  
131 137 //shift and add the new value to the array
132   - j[0] = j[1];
133   - j[1] = jnew;
  138 + b[0] = b[1];
  139 + b[1] = bnew;
134 140 }
135 141  
136 142  
... ...