00001 // This is mul/vpdfl/vpdfl_gaussian_kernel_pdf_sampler.cxx 00002 #ifdef VCL_NEEDS_PRAGMA_INTERFACE 00003 #pragma implementation 00004 #endif 00005 //: 00006 // \file 00007 // \author Tim Cootes 00008 // \brief Sampler class for gaussian kernel PDF. 00009 00010 #include "vpdfl_gaussian_kernel_pdf_sampler.h" 00011 00012 #include <vcl_cassert.h> 00013 #include <vpdfl/vpdfl_gaussian_kernel_pdf.h> 00014 00015 //======================================================================= 00016 // Dflt ctor 00017 //======================================================================= 00018 00019 vpdfl_gaussian_kernel_pdf_sampler::vpdfl_gaussian_kernel_pdf_sampler(): 00020 rng_(9667566ul) 00021 { 00022 } 00023 00024 //======================================================================= 00025 // Destructor 00026 //======================================================================= 00027 00028 vpdfl_gaussian_kernel_pdf_sampler::~vpdfl_gaussian_kernel_pdf_sampler() 00029 { 00030 } 00031 00032 00033 //======================================================================= 00034 // Method: is_a 00035 //======================================================================= 00036 00037 vcl_string vpdfl_gaussian_kernel_pdf_sampler::is_a() const 00038 { 00039 return vcl_string("vpdfl_gaussian_kernel_pdf_sampler"); 00040 } 00041 00042 //======================================================================= 00043 // Method: is_class 00044 //======================================================================= 00045 00046 bool vpdfl_gaussian_kernel_pdf_sampler::is_class(vcl_string const& s) const 00047 { 00048 return vpdfl_sampler_base::is_class(s) || s==vpdfl_gaussian_kernel_pdf_sampler::is_a(); 00049 } 00050 00051 //======================================================================= 00052 // Method: clone 00053 //======================================================================= 00054 00055 vpdfl_sampler_base* vpdfl_gaussian_kernel_pdf_sampler::clone() const 00056 { 00057 return new vpdfl_gaussian_kernel_pdf_sampler(*this); 00058 } 00059 //======================================================================= 00060 00061 void vpdfl_gaussian_kernel_pdf_sampler::reseed(unsigned long seed) 00062 { 00063 rng_.reseed(seed); 00064 } 00065 //======================================================================= 00066 00067 //: Set model for which this is an instance 00068 // Error check that it is an axis gaussian. 00069 void vpdfl_gaussian_kernel_pdf_sampler::set_model(const vpdfl_pdf_base& model) 00070 { 00071 assert(model.is_class("vpdfl_gaussian_kernel_pdf")); 00072 // cannot use dynamic_cast<> without rtti - PVr 00073 // rtti currently turned off 00074 vpdfl_sampler_base::set_model(model); 00075 } 00076 00077 void vpdfl_gaussian_kernel_pdf_sampler::sample_component(vnl_vector<double>& x, 00078 int j) 00079 { 00080 const vpdfl_gaussian_kernel_pdf & kpdf = static_cast<const vpdfl_gaussian_kernel_pdf &>(model()); 00081 00082 int n_dims = kpdf.n_dims(); 00083 x.set_size(n_dims); 00084 00085 const double* m = kpdf.centre()[j].data_block(); 00086 double w = kpdf.width()[j]; 00087 00088 double* x_data = x.data_block(); 00089 for (int i=0;i<n_dims;++i) 00090 x_data[i] = m[i] + w*rng_.normal(); 00091 } 00092 00093 00094 //======================================================================= 00095 00096 void vpdfl_gaussian_kernel_pdf_sampler::sample(vnl_vector<double>& x) 00097 { 00098 const vpdfl_gaussian_kernel_pdf & kpdf = static_cast<const vpdfl_gaussian_kernel_pdf &>(model()); 00099 int n = kpdf.centre().size(); 00100 00101 // Select component 00102 int j = 0; 00103 if (n>1) j=rng_.lrand32(0,n-1); 00104 00105 sample_component(x,j); 00106 } 00107 00108 //: Fill x with samples possibly chosen so as to represent the distribution 00109 // Sample sequentially from each component. 00110 void vpdfl_gaussian_kernel_pdf_sampler::regular_samples( 00111 vcl_vector<vnl_vector<double> >& x) 00112 { 00113 const vpdfl_gaussian_kernel_pdf & kpdf = static_cast<const vpdfl_gaussian_kernel_pdf &>(model()); 00114 int n_k = kpdf.centre().size(); 00115 00116 int n_samples = x.size(); 00117 00118 for (int i=0;i<n_samples;++i) 00119 { 00120 // Ensure sample at centre. (This biases towards the centres somewhat) 00121 // if (i<n_k) 00122 // x[i]=kpdf.centre()[i]; 00123 // else 00124 sample_component(x[i],i%n_k); 00125 } 00126 } 00127 00128 //======================================================================= 00129 00130 00131 //: Return a reference to the pdf model 00132 // This is properly cast. 00133 const vpdfl_gaussian_kernel_pdf& vpdfl_gaussian_kernel_pdf_sampler::gaussian_kernel_pdf() const 00134 { 00135 return static_cast<const vpdfl_gaussian_kernel_pdf&>( model()); 00136 } 00137