00001 #ifndef CCTBX_MAPTBX_PEAK_SEARCH_H
00002 #define CCTBX_MAPTBX_PEAK_SEARCH_H
00003
00004 #include <scitbx/indexed_value.h>
00005 #include <scitbx/array_family/selections.h>
00006 #include <scitbx/array_family/accessors/c_grid.h>
00007 #include <scitbx/array_family/accessors/c_grid_padded.h>
00008 #include <scitbx/array_family/loops.h>
00009 #include <scitbx/array_family/sort.h>
00010 #include <scitbx/array_family/tiny_algebra.h>
00011 #include <scitbx/sym_mat3.h>
00012 #include <cctbx/error.h>
00013 #include <scitbx/math/modulo.h>
00014 #include <cctbx/import_scitbx_af.h>
00015
00016 namespace cctbx { namespace maptbx {
00017
00019
00041 template <typename DataType,
00042 typename TagType>
00043 void
00044 peak_search_unit_cell(
00045 af::const_ref<DataType, af::c_grid_padded<3> > const& data,
00046 af::ref<TagType, af::c_grid<3> > const& tags,
00047 int level)
00048 {
00049 CCTBX_ASSERT(tags.accessor().all_eq(data.accessor().focus()));
00050 CCTBX_ASSERT(!data.accessor().is_padded());
00051 const DataType* pdata = data.begin();
00052 TagType* ptags = tags.begin();
00053 int nx = tags.accessor()[0];
00054 int ny = tags.accessor()[1];
00055 int nz = tags.accessor()[2];
00056 int nk = nz;
00057 int nj_nk = ny * nz;
00058 int ni_nj_nk = nx * ny * nz;
00059
00060 int im, jm, km;
00061 int i0, j0, k0;
00062 int ip, jp, kp;
00063 int ibreak, jbreak, kbreak;
00064
00065
00066 for (int i = 0; i < ni_nj_nk; i++) {
00067 if (ptags[i] < 0) ptags[i] = -1;
00068 }
00069
00070 const DataType* pd = pdata;
00071 TagType* pf = ptags;
00072 im = ni_nj_nk - nj_nk;
00073 i0 = 0;
00074 ip = nj_nk;
00075 ibreak = ni_nj_nk;
00076 while (ip < ibreak) {
00077 jm = nj_nk - nk;
00078 j0 = 0;
00079 jp = nk;
00080 jbreak = nj_nk;
00081 while (jp < jbreak) {
00082 km = nk - 1;
00083 k0 = 0;
00084 kp = 1;
00085 kbreak = nk;
00086 while (kp < kbreak) {
00087 TagType* indep_pf;
00088 if (*pf < 0) indep_pf = pf;
00089 else indep_pf = &ptags[*pf];
00090 while (! (*indep_pf < -1)) {
00091 if (level >= 1) {
00092
00093
00094
00095 if (*pd < pdata[im + j0 + k0]) break;
00096 if (*pd < pdata[ip + j0 + k0]) break;
00097 if (*pd < pdata[i0 + jm + k0]) break;
00098 if (*pd < pdata[i0 + jp + k0]) break;
00099 if (*pd < pdata[i0 + j0 + km]) break;
00100 if (*pd < pdata[i0 + j0 + kp]) break;
00101 if (level >= 2) {
00102
00103
00104
00105 if (*pd < pdata[im + jm + k0]) break;
00106 if (*pd < pdata[ip + jp + k0]) break;
00107 if (*pd < pdata[im + j0 + km]) break;
00108 if (*pd < pdata[ip + j0 + kp]) break;
00109 if (*pd < pdata[i0 + jm + km]) break;
00110 if (*pd < pdata[i0 + jp + kp]) break;
00111 if (*pd < pdata[im + jp + k0]) break;
00112 if (*pd < pdata[ip + jm + k0]) break;
00113 if (*pd < pdata[im + j0 + kp]) break;
00114 if (*pd < pdata[ip + j0 + km]) break;
00115 if (*pd < pdata[i0 + jm + kp]) break;
00116 if (*pd < pdata[i0 + jp + km]) break;
00117 if (level >= 3) {
00118
00119
00120
00121 if (*pd < pdata[im + jm + km]) break;
00122 if (*pd < pdata[ip + jp + kp]) break;
00123 if (*pd < pdata[im + jm + kp]) break;
00124 if (*pd < pdata[ip + jp + km]) break;
00125 if (*pd < pdata[im + jp + km]) break;
00126 if (*pd < pdata[ip + jm + kp]) break;
00127 if (*pd < pdata[im + jp + kp]) break;
00128 if (*pd < pdata[ip + jm + km]) break;
00129 }
00130 }
00131 }
00132 *indep_pf = -2;
00133 break;
00134 }
00135 pd++;
00136 pf++;
00137 km = k0;
00138 k0 = kp;
00139 kp++;
00140 if (kp == nk) { kp = 0; kbreak = 1; }
00141 }
00142 jm = j0;
00143 j0 = jp;
00144 jp += nk;
00145 if (jp == nj_nk) { jp = 0; jbreak = nk; }
00146 }
00147 im = i0;
00148 i0 = ip;
00149 ip += nj_nk;
00150 if (ip == ni_nj_nk) { ip = 0; ibreak = nj_nk; }
00151 }
00152 }
00153
00154 template <typename ValueType,
00155 typename CountType = long>
00156 struct peak_histogram
00157 {
00158 peak_histogram() {}
00159
00160 template <typename DataType,
00161 typename TagType>
00162 peak_histogram(
00163 af::const_ref<DataType, af::c_grid_padded<3> > const& data,
00164 af::ref<TagType, af::c_grid<3> > const& tags,
00165 std::size_t n_slots=1000)
00166 :
00167 slots(n_slots)
00168 {
00169 CCTBX_ASSERT(data.accessor().focus().all_eq(tags.accessor()));
00170 CCTBX_ASSERT(n_slots > 0);
00171 std::size_t i;
00172 for(i=0;i<data.size();i++) {
00173 if (tags[i] == -2) break;
00174 }
00175 if (i == data.size()) {
00176 data_min = 0;
00177 data_max = 0;
00178 }
00179 else {
00180 data_min = data[i];
00181 data_max = data[i];
00182 for(i++;i<data.size();i++) {
00183 if (tags[i] != -2) continue;
00184 if (data_min > data[i]) data_min = data[i];
00185 if (data_max < data[i]) data_max = data[i];
00186 }
00187 }
00188 slot_width = (data_max - data_min) / slots.size();
00189 slots.fill(0);
00190 for(i=0;i<data.size();i++) {
00191 if (tags[i] != -2) continue;
00192 ValueType d = data[i] - data_min;
00193 std::size_t i_slot = 0;
00194 if (d != 0 && d >= slot_width) {
00195 i_slot = std::size_t(d / slot_width);
00196 if (i_slot >= slots.size()) i_slot = slots.size() - 1;
00197 }
00198 slots[i_slot]++;
00199 }
00200 }
00201
00202 ValueType
00203 get_cutoff(CountType const& max_peaks,
00204 ValueType const& tolerance=1.e-4) const
00205 {
00206 CountType cum = 0;
00207 std::size_t i = slots.size();
00208 for (; i; i--) {
00209 cum += slots[i-1];
00210 if (cum > max_peaks) break;
00211 }
00212 return data_min + i * slot_width + slot_width * tolerance;
00213 }
00214
00215 ValueType data_min;
00216 ValueType data_max;
00217 ValueType slot_width;
00218 af::shared<CountType> slots;
00219 };
00220
00221 template <typename GridIndexType = af::tiny<long, 3>,
00222 typename SiteType = scitbx::vec3<double>,
00223 typename ValueType = double>
00224 class peak_list
00225 {
00226 public:
00227 peak_list() {}
00228
00229 template <typename DataType,
00230 typename TagType>
00231 peak_list(
00232 af::const_ref<DataType, af::c_grid_padded<3> > const& data,
00233 af::ref<TagType, af::c_grid<3> > const& tags,
00234 int peak_search_level=1,
00235 std::size_t max_peaks=0,
00236 bool interpolate=true)
00237 :
00238 gridding_(data.accessor().focus())
00239 {
00240 peak_search_unit_cell(data, tags, peak_search_level);
00241 ValueType peak_cutoff = 0;
00242 bool use_cutoff = false;
00243 if (max_peaks) {
00244 peak_histogram<ValueType> hist(data, tags);
00245 peak_cutoff = hist.get_cutoff(max_peaks);
00246 use_cutoff = true;
00247 }
00248 process_peaks(data, tags, peak_cutoff, use_cutoff, interpolate);
00249 }
00250
00251 template <typename DataType,
00252 typename TagType>
00253 peak_list(
00254 af::const_ref<DataType, af::c_grid_padded<3> > const& data,
00255 af::ref<TagType, af::c_grid<3> > const& tags,
00256 int peak_search_level,
00257 ValueType peak_cutoff,
00258 std::size_t max_peaks,
00259 bool interpolate=true)
00260 :
00261 gridding_(data.accessor().focus())
00262 {
00263 peak_search_unit_cell(data, tags, peak_search_level);
00264 if (max_peaks) {
00265 peak_histogram<ValueType> hist(data, tags);
00266 peak_cutoff = std::max(peak_cutoff, hist.get_cutoff(max_peaks));
00267 }
00268 process_peaks(data, tags, peak_cutoff, true, interpolate);
00269 }
00270
00271 GridIndexType const&
00272 gridding() const { return gridding_; }
00273
00274 std::size_t
00275 size() const
00276 {
00277 CCTBX_ASSERT(grid_heights().size() == grid_indices().size());
00278 CCTBX_ASSERT(sites().size() == grid_indices().size());
00279 CCTBX_ASSERT(heights().size() == grid_indices().size());
00280 return grid_indices().size();
00281 }
00282
00285 af::shared<GridIndexType>
00286 grid_indices() const { return grid_indices_; }
00287
00288 af::shared<ValueType>
00289 grid_heights() const { return grid_heights_; }
00290
00291 af::shared<SiteType>
00292 sites() const { return sites_; }
00293
00294 af::shared<ValueType>
00295 heights() const { return heights_; }
00296
00297 protected:
00298 GridIndexType gridding_;
00299 af::shared<GridIndexType> grid_indices_;
00300 af::shared<ValueType> grid_heights_;
00301 af::shared<SiteType> sites_;
00302 af::shared<ValueType> heights_;
00303
00304 template <typename DataType,
00305 typename TagType>
00306 void
00307 process_peaks(
00308 af::const_ref<DataType, af::c_grid_padded<3> > const& data,
00309 af::const_ref<TagType, af::c_grid<3> > const& tags,
00310 ValueType const& cutoff,
00311 bool use_cutoff,
00312 bool interpolate)
00313 {
00314 collect_peaks(data, tags, cutoff, use_cutoff);
00315 if (interpolate) interpolate_sites_and_heights(data);
00316 else copy_sites_and_heights();
00317 sort();
00318 }
00319
00320 template <typename DataType,
00321 typename TagType>
00322 void
00323 collect_peaks(
00324 af::const_ref<DataType, af::c_grid_padded<3> > const& data,
00325 af::const_ref<TagType, af::c_grid<3> > const& tags,
00326 ValueType const& cutoff,
00327 bool use_cutoff)
00328 {
00329 typedef typename af::c_grid<3>::index_type index_type;
00330 af::nested_loop<index_type> loop(data.accessor().focus());
00331 for (index_type const& pivot = loop(); !loop.over(); loop.incr()) {
00332 if (tags(pivot) == -2 && (!use_cutoff || data(pivot) >= cutoff)) {
00333 grid_indices_.push_back(pivot);
00334 grid_heights_.push_back(data(pivot));
00335 }
00336 }
00337 }
00338
00339 void
00340 copy_sites_and_heights()
00341 {
00342 af::const_ref<GridIndexType> gi = grid_indices_.const_ref();
00343 af::tiny<ValueType, 3> gr(gridding_);
00344 sites_.reserve(gi.size());
00345 for(std::size_t i=0;i<gi.size();i++) {
00346 sites_.push_back(af::tiny<ValueType, 3>(gi[i]) / gr);
00347 }
00348 heights_.assign(grid_heights_);
00349 }
00350
00351 template <typename DataType>
00352 void
00353 interpolate_sites_and_heights(
00354 af::const_ref<DataType, af::c_grid_padded<3> > const& data,
00355 ValueType epsilon=1.e-6)
00356 {
00357 af::const_ref<GridIndexType> gi = grid_indices_.const_ref();
00358 af::const_ref<ValueType> gh = grid_heights_.const_ref();
00359 af::tiny<ValueType, 3> gr(gridding_);
00360 sites_.reserve(gi.size());
00361 heights_.reserve(gi.size());
00362 for(std::size_t i_peak=0;i_peak<gi.size();i_peak++) {
00363 af::tiny<ValueType, 3> site = gi[i_peak];
00364 ValueType height = gh[i_peak];
00365 long u0 = gi[i_peak][0];
00366 long up = scitbx::math::mod_positive(u0+1, gridding_[0]);
00367 long um = scitbx::math::mod_positive(u0-1, gridding_[0]);
00368 long v0 = gi[i_peak][1];
00369 long vp = scitbx::math::mod_positive(v0+1, gridding_[1]);
00370 long vm = scitbx::math::mod_positive(v0-1, gridding_[1]);
00371 long w0 = gi[i_peak][2];
00372 long wp = scitbx::math::mod_positive(w0+1, gridding_[2]);
00373 long wm = scitbx::math::mod_positive(w0-1, gridding_[2]);
00374 scitbx::vec3<ValueType> b(
00375 (data(um,v0,w0) - data(up,v0,w0)) / 2,
00376 (data(u0,vm,w0) - data(u0,vp,w0)) / 2,
00377 (data(u0,v0,wm) - data(u0,v0,wp)) / 2);
00378 ValueType two_gh = 2 * height;
00379 scitbx::sym_mat3<ValueType> a(
00380 data(um,v0,w0) + data(up,v0,w0) - two_gh,
00381 data(u0,vm,w0) + data(u0,vp,w0) - two_gh,
00382 data(u0,v0,wm) + data(u0,v0,wp) - two_gh,
00383 ( (data(up,vp,w0) + data(um,vm,w0))
00384 - (data(up,vm,w0) + data(um,vp,w0))) / 4,
00385 ( (data(up,v0,wp) + data(um,v0,wm))
00386 - (data(up,v0,wm) + data(um,v0,wp))) / 4,
00387 ( (data(u0,vp,wp) + data(u0,vm,wm))
00388 - (data(u0,vp,wm) + data(u0,vm,wp))) / 4);
00389 scitbx::sym_mat3<ValueType> a_inv = a.co_factor_matrix_transposed();
00390 ValueType a_inv_max = af::max_absolute(a_inv);
00391 ValueType a_det = a.determinant();
00392 if (scitbx::fn::absolute(a_det) > a_inv_max * epsilon) {
00393 a_inv /= a_det;
00394 af::tiny<ValueType, 3> x = a_inv * b;
00395 if (af::max_absolute(x) < 1) {
00396 site += x;
00397 height -= b * scitbx::vec3<ValueType>(x);
00398 for(std::size_t i=0;i<3;i++) {
00399 height += a[i]*x[i]*x[i]/2;
00400 }
00401 height += a[3]*x[0]*x[1]
00402 + a[4]*x[0]*x[2]
00403 + a[5]*x[1]*x[2];
00404 }
00405 }
00406 site /= gr;
00407 sites_.push_back(site);
00408 heights_.push_back(height);
00409 }
00410 }
00411
00412 void
00413 sort()
00414 {
00415 af::shared<std::size_t>
00416 perm = af::sort_permutation(heights_.const_ref(), true);
00417 af::const_ref<std::size_t> p = perm.const_ref();
00418 grid_indices_ = af::select(grid_indices_.const_ref(), p);
00419 grid_heights_ = af::select(grid_heights_.const_ref(), p);
00420 sites_ = af::select(sites_.const_ref(), p);
00421 heights_ = af::select(heights_.const_ref(),p);
00422 }
00423 };
00424
00425 }}
00426
00427 #endif // CCTBX_MAPTBX_PEAK_SEARCH_H
