<commit_before>/*
 * image_processing.h
 *
 *  Created on: Feb 28, 2010
 *      Author: <NAME>
 */

#ifndef IMAGE_PROCESSING_H_
#define IMAGE_PROCESSING_H_
 
#include "Tree.h"
#ifdef ENABLE_FFTW
#include "fftw3.h"
#endif


// forward declaration
struct Grid;

/** \ingroup Image
 * \brief Takes image structure and pixelizes the flux into regular pixel grid which then
 * can be exported as a fits file, smoothed, etc. like an image.
 *
 */
class PixelMap
{
public:
	PixelMap(const PixelMap& pmap, double res_ratio);
	PixelMap();
	PixelMap(const PixelMap& other);
	PixelMap(const PixelMap& pmap, const double* center, std::size_t Npixels);
	PixelMap(const double* center, std::size_t Npixels, double resolution);
	PixelMap(const double* center, std::size_t Nx, std::size_t Ny, double resolution);
	PixelMap(std::string fitsfilename,double resolution = -1);
	~PixelMap();
	
	PixelMap& operator=(PixelMap other);
	
	inline bool valid() const { return map.size(); };
	inline std::size_t size() const { return map.size(); };
	
	inline std::size_t getNx() const { return Nx; }
	inline std::size_t getNy() const { return Ny; }
	inline double getRangeX() const { return rangeX; }
	inline double getRangeY() const { return rangeY; }
	inline const double* getCenter() const{ return center; }
	inline double getResolution() const { return resolution; }
	
	void Clean();

	void AddImages(ImageInfo *imageinfo,int Nimages,float rescale = 1.);
	void AddCurve(ImageInfo *curve,double value);
	void drawline(double x1[],double x2[],double value);
  void drawcircle(PosType r_center[],PosType radius,PosType value);
	void AddGrid(Grid &grid,double value = 1.0);
  void AddGrid(Grid &grid,LensingVariable val);

	void Renormalize(double factor);
	void AddValue(std::size_t i, double value);
	void AssignValue(std::size_t i, double value);
	void printASCII() const;
	void printASCIItoFile(std::string filename) const;
	void printFITS(std::string filename, bool verbose = false) const;

	void smooth(double sigma);

	inline double getValue(std::size_t i) const { return map[i]; }
	inline double & operator[](std::size_t i) { return map[i]; };
	inline double operator()(std::size_t i) const { return map[i]; };
	
	PixelMap& operator+=(const PixelMap& rhs);
	friend PixelMap operator+(const PixelMap&, const PixelMap&);

	PixelMap& operator-=(const PixelMap& rhs);
	friend PixelMap operator-(const PixelMap&, const PixelMap&);
    
	PixelMap& operator*=(const PixelMap& rhs);
	friend PixelMap operator*(const PixelMap&, const PixelMap&);

	PixelMap& operator*=(PosType b);
	friend PixelMap operator*(const PixelMap&, PosType b);
	
	const std::valarray<double>& data() const { return map; }
	
	bool agrees(const PixelMap& other) const;
	
	friend void swap(PixelMap&, PixelMap&);
  
  /// return average pixel value
  PosType ave() const;
  /// Total number of pixels
  size_t size(){return map.size();}
	
  void FindArc(PosType &radius,PosType *xc,PosType *arc_center,PosType &arclength,PosType &width
                         ,PosType threshold);
  
  /// get the index for a position, returns -1 if out of map, this version returns the 2D grid coordinates
  long find_index(PosType const x[],long &ix,long &iy);
  /// get the index for a position, returns -1 if out of map
  long find_index(PosType const x[]);
  
  /// get the index for a position, returns -1 if out of map, this version returns the 2D grid coordinates
  long find_index(PosType const x,PosType const y,long &ix,long &iy);
  /// get the index for a position, returns -1 if out of map
  long find_index(PosType const x,PosType const y);
  
  /// get the index for a position, returns -1 if out of map
  void find_position(PosType x[],std::size_t const index);
  /// get the index for a position, returns -1 if out of map
  void find_position(PosType x[],std::size_t const ix,std::size_t const iy);
  
  /// interpolate to point x[]
  PosType linear_interpolate(PosType x[]);
  
#ifdef ENABLE_FFTW

  void PowerSpectrum(std::vector<PosType> &power_spectrum,bool logbins){
    
    int Nbins = power_spectrum.size();
    std::vector<size_t> NinBins(Nbins,0);
    std::vector<size_t> bin_vec(Nbins + 1);
    
    PosType nyquist = sqrt(Nx*Nx + Ny*Ny)/2;
    
    int l;
    if(!logbins) for(l=0;l<Nbins+1;++l) bin_vec[l] = l*nyquist/Nbins;
    else{
      PosType f = pow(nyquist,1.0/Nbins);
      bin_vec[0]=0;
      bin_vec[1]=f;
      for(l=2;l<Nbins+1;++l) bin_vec[l] = bin_vec[l-1]*f;
    }
    
    //if(Nx != Ny ) throw std::runtime_error("PixelMap must be square to find Power Spectrum");
    
    //PosType *tmp_map = new PosType[map.size()];
    //for(size_t i=0; i<map.size();++i ) tmp_map[i] = map[i];
    
    // Fourier tranform the plane
    fftw_complex *fftoutput = new fftw_complex[Nx*(Ny/2+1)];
    //fftw_plan p = fftw_plan_dft_r2c_2d(Nx,Ny,tmp_map,fftoutput,FFTW_ESTIMATE);
    fftw_plan p = fftw_plan_dft_r2c_2d(Nx,Ny,std::begin(map),fftoutput,FFTW_ESTIMATE);
    fftw_execute( p );
    fftw_destroy_plan(p);
    //delete[] tmp_map;

    for(l=0;l<Nbins;++l) power_spectrum[l] = 0.0;

    // build modes for each pixel in the fourier space
    size_t k2,ii;
    long kx;
    for( size_t i=0; i<Nx; i++ ){
      kx=(i<Nx/2)? (i):(i-Nx);
      for( int j=0; j<Ny/2+1; j++ ){
        k2 = sqrt(kx*kx + j*j);
        ii = j+(Ny/2+1)*i;
        
        l = Utilities::locate<size_t>(bin_vec, k2);
        //std::cout << fftoutput[ii][0] << "  " << fftoutput[ii][1] << std::endl;
        
        if(l >= 0 && l < Nbins ) power_spectrum[l] += fftoutput[ii][0]*fftoutput[ii][0] +
                                                     fftoutput[ii][1]*fftoutput[ii][1];
        ++(NinBins[l]);
      }
    }
    for(l=0;l<Nbins;++l) power_spectrum[l] /= NinBins[l];
    
    delete[] fftoutput;
  }
#endif

private:
	std::valarray<double> map;
  void AddGrid_(PointList list,LensingVariable val);

	std::size_t Nx;
	std::size_t Ny;
	double resolution,rangeX,rangeY,center[2];
	double map_boundary_p1[2],map_boundary_p2[2];

	double LeafPixelArea(IndexType i,Branch * branch1);
	void PointsWithinLeaf(Branch * branch1, std::list <unsigned long> &neighborlist);
	bool inMapBox(Branch * branch1) const;
	bool inMapBox(double * branch1) const;
  
};
<commit_msg>Work on PixelMap::PowerSpectrum()<commit_after>/*
 * image_processing.h
 *
 *  Created on: Feb 28, 2010
 *      Author: <NAME>
 */

#ifndef IMAGE_PROCESSING_H_
#define IMAGE_PROCESSING_H_
 
#include "Tree.h"
#ifdef ENABLE_FFTW
#include "fftw3.h"
#endif


// forward declaration
struct Grid;

/** \ingroup Image
 * \brief Takes image structure and pixelizes the flux into regular pixel grid which then
 * can be exported as a fits file, smoothed, etc. like an image.
 *
 */
class PixelMap
{
public:
	PixelMap(const PixelMap& pmap, double res_ratio);
	PixelMap();
	PixelMap(const PixelMap& other);
	PixelMap(const PixelMap& pmap, const double* center, std::size_t Npixels);
	PixelMap(const double* center, std::size_t Npixels, double resolution);
	PixelMap(const double* center, std::size_t Nx, std::size_t Ny, double resolution);
	PixelMap(std::string fitsfilename,double resolution = -1);
	~PixelMap();
	
	PixelMap& operator=(PixelMap other);
	
	inline bool valid() const { return map.size(); };
	inline std::size_t size() const { return map.size(); };
	
	inline std::size_t getNx() const { return Nx; }
	inline std::size_t getNy() const { return Ny; }
	inline double getRangeX() const { return rangeX; }
	inline double getRangeY() const { return rangeY; }
	inline const double* getCenter() const{ return center; }
	inline double getResolution() const { return resolution; }
	
	void Clean();

	void AddImages(ImageInfo *imageinfo,int Nimages,float rescale = 1.);
	void AddCurve(ImageInfo *curve,double value);
	void drawline(double x1[],double x2[],double value);
  void drawcircle(PosType r_center[],PosType radius,PosType value);
	void AddGrid(Grid &grid,double value = 1.0);
  void AddGrid(Grid &grid,LensingVariable val);

	void Renormalize(double factor);
	void AddValue(std::size_t i, double value);
	void AssignValue(std::size_t i, double value);
	void printASCII() const;
	void printASCIItoFile(std::string filename) const;
	void printFITS(std::string filename, bool verbose = false) const;

	void smooth(double sigma);

	inline double getValue(std::size_t i) const { return map[i]; }
	inline double & operator[](std::size_t i) { return map[i]; };
	inline double operator()(std::size_t i) const { return map[i]; };
	
	PixelMap& operator+=(const PixelMap& rhs);
	friend PixelMap operator+(const PixelMap&, const PixelMap&);

	PixelMap& operator-=(const PixelMap& rhs);
	friend PixelMap operator-(const PixelMap&, const PixelMap&);
    
	PixelMap& operator*=(const PixelMap& rhs);
	friend PixelMap operator*(const PixelMap&, const PixelMap&);

	PixelMap& operator*=(PosType b);
	friend PixelMap operator*(const PixelMap&, PosType b);
	
	const std::valarray<double>& data() const { return map; }
	
	bool agrees(const PixelMap& other) const;
	
	friend void swap(PixelMap&, PixelMap&);
  
  /// return average pixel value
  PosType ave() const;
  /// Total number of pixels
  size_t size(){return map.size();}
	
  void FindArc(PosType &radius,PosType *xc,PosType *arc_center,PosType &arclength,PosType &width
                         ,PosType threshold);
  
  /// get the index for a position, returns -1 if out of map, this version returns the 2D grid coordinates
  long find_index(PosType const x[],long &ix,long &iy);
  /// get the index for a position, returns -1 if out of map
  long find_index(PosType const x[]);
  
  /// get the index for a position, returns -1 if out of map, this version returns the 2D grid coordinates
  long find_index(PosType const x,PosType const y,long &ix,long &iy);
  /// get the index for a position, returns -1 if out of map
  long find_index(PosType const x,PosType const y);
  
  /// get the index for a position, returns -1 if out of map
  void find_position(PosType x[],std::size_t const index);
  /// get the index for a position, returns -1 if out of map
  void find_position(PosType x[],std::size_t const ix,std::size_t const iy);
  
  /// interpolate to point x[]
  PosType linear_interpolate(PosType x[]);
  
#ifdef ENABLE_FFTW

  void PowerSpectrum(std::vector<PosType> &power_spectrum,std::vector<PosType> &lvec,bool logbins){
    
    int Nbins = power_spectrum.size();
    std::vector<size_t> NinBins(Nbins,0);
    std::vector<size_t> bin_vec(Nbins + 1);
    
    PosType nyquist = sqrt(Nx*Nx + Ny*Ny)/2;
    
    int l;
    if(!logbins) for(l=0;l<Nbins+1;++l) bin_vec[l] = l*nyquist/Nbins;
    else{
      PosType f = pow(nyquist,1.0/Nbins);
      bin_vec[0]=0;
      bin_vec[1]=f;
      for(l=2;l<Nbins+1;++l) bin_vec[l] = bin_vec[l-1]*f;
    }
    
    //if(Nx != Ny ) throw std::runtime_error("PixelMap must be square to find Power Spectrum");
    
    //PosType *tmp_map = new PosType[map.size()];
    //for(size_t i=0; i<map.size();++i ) tmp_map[i] = map[i];
    
    // Fourier tranform the plane
    fftw_complex *fftoutput = new fftw_complex[Nx*(Ny/2+1)];
    //fftw_plan p = fftw_plan_dft_r2c_2d(Nx,Ny,tmp_map,fftoutput,FFTW_ESTIMATE);
    fftw_plan p = fftw_plan_dft_r2c_2d(Nx,Ny,std::begin(map),fftoutput,FFTW_ESTIMATE);
    fftw_execute( p );
    fftw_destroy_plan(p);
    //delete[] tmp_map;

    for(l=0;l<Nbins;++l) power_spectrum[l] = 0.0;

    // build modes for each pixel in the fourier space
    size_t k2,ii;
    long kx;
    for( size_t i=0; i<Nx; i++ ){
      kx=(i<Nx/2)? (i):(i-Nx);
      for( int j=0; j<Ny/2+1; j++ ){
        k2 = sqrt(kx*kx + j*j);
        ii = j+(Ny/2+1)*i;
        
        l = Utilities::locate<size_t>(bin_vec, k2);
        //std::cout << fftoutput[ii][0] << "  " << fftoutput[ii][1] << std::endl;
        
        if(l >= 0 && l < Nbins ) power_spectrum[l] += fftoutput[ii][0]*fftoutput[ii][0] +
                                                     fftoutput[ii][1]*fftoutput[ii][1];
        ++(NinBins[l]);
      }
    }
    for(l=0;l<Nbins;++l) power_spectrum[l] /= NinBins[l];
    
    lvec.resize(Nbins);
    if(!logbins) for(l=0;l<Nbins;++l) lvec[l] = (bin_vec[l] + bin_vec[l+1])/2/resolution/nyquist;
    
    delete[] fftoutput;
  }
#endif

private:
	std::valarray<double> map;
  void AddGrid_(PointList list,LensingVariable val);

	std::size_t Nx;
	std::size_t Ny;
	double resolution,rangeX,rangeY,center[2];
	double map_boundary_p1[2],map_boundary_p2[2];

	double LeafPixelArea(IndexType i,Branch * branch1);
	void PointsWithinLeaf(Branch * branch1, std::list <unsigned long> &neighborlist);
	bool inMapBox(Branch * branch1) const;
	bool inMapBox(double * branch1) const;
  
};

{'fraction_non_alphanumeric': 0.138664128192924, 'fraction_numerical': 0.011502459146438203, 'mean_word_length': 3.8845408756296007, 'pattern_counts': {'":': 0, '<': 55, '<?xml version=': 0, '>': 22, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '89932', 'n_tokens_mistral': 4676, 'n_tokens_neox': 4470, 'n_words': 1193}