source: sans/Dev/trunk/NCNR_User_Procedures/Reduction/VSANS/V_TubeAdjustments.ipf @ 965

Last change on this file since 965 was 965, checked in by srkline, 7 years ago

intermediate changes to the VSANS r/w routines, a utility for processing ISO8601 dates, and some testing procedures for working with the banks of tubes.

File size: 13.4 KB
Line 
1#pragma rtGlobals=3             // Use modern global access method and strict wave access.
2
3
4
5
6
7// the main routines are:
8
9//(1)
10//to get from individual tubes to an array
11//      Tubes_to_Array()                       
12
13//(2)
14// then to locate all of the peak positions
15//      MakeTableForPeaks(numTube,numPeak)             
16//      Identify_AllPeaks()
17//              AutoFindPeaksCustom()           // if Identify_AllPeaks  doesn't work -try this, setting the "noise" to 1 and smoothing to 2
18
19//(3)
20// fit to find all of the quadratic coefficients
21//      MakeTableForFitCoefs(numTube,numCoef)
22//      PlotFit_AllPeaks()
23
24
25//(4)
26// then pick a display method
27//
28//      Make_mm_tubes()
29//      Append_Adjusted_mm()
30//
31//      MakeMatrix_PixCenter()
32//      FillMatrix_Pix_Center(ind)
33//
34//
35// -or- (note that the pack_image wave that is generated here is for display ONLY)
36// --since it is interpolated
37//
38// Interpolate_mm_tubes()
39
40
41// The function most used externally is:
42// V_TubePix_to_mm(c0,c1,c2,pix)
43//
44// which will return the real space distance (in mm?) for a given pixel
45// and the tube's coefficients. The distance is relative to the zero position on the
46// detector (which was defined when the coefficients were determined)
47
48
49
50
51
52
53
54// (0) -- what I start with:
55// -- a table of the mm spacing of the slots (20 of them)
56// -- masked data from each of the (8) tubes
57// -- the table of slots may need to be corrected for parallax, depending on the geometry of the test
58// ** In the table of slots, pick a slot near the center, and SET that to ZERO. Then all of the other
59//   distances are relative to that zero point. This is a necessary reference point.
60//
61
62
63// (1) -- get the individual tubes into an array
64//
65//
66Proc Tubes_to_Array()
67        Make/O/D/N=(8,1127) pack
68        edit pack
69        display;appendimage pack
70        pack[0][] = tube1[q]
71        pack[1][] = tube2[q]
72        pack[2][] = tube3[q]
73        pack[3][] = tube4[q]
74        pack[4][] = tube5[q]
75        pack[5][] = tube6[q]
76        pack[6][] = tube7[q]
77        pack[7][] = tube8[q]
78        ModifyImage pack ctab= {*,*,ColdWarm,0}
79End
80
81// (2) -- for each of the tubes, find the x-position (in pixels) of each of the (20) peaks
82// -- load the Analysis Package "MultiPeakFit 2"
83//
84// automatically find the peaks (after including MultiPeakFit 2)
85//              AutomaticallyFindPeaks()
86//
87//-- or if having difficulty
88//              AutoFindPeaksCustom()           // try this, setting the "noise" to 1 and smoothing to 2
89//
90// -- if really having difficulty, you can do the "full" MultiPeak Fit
91//
92// -- If (hopefully) using the easy way, the results are in:
93// root:WA_PeakCentersY,root:WA_PeakCentersX
94//
95// -- so to see the results:
96//¥Edit/K=0  root:WA_PeakCentersY,root:WA_PeakCentersX
97//
98// -- then sort the results - they seem to be in no real order...
99//¥Sort WA_PeakCentersX WA_PeakCentersY,WA_PeakCentersX
100//
101Proc MakeTableForPeaks(numTube,numPeak)
102        Variable numTube,numPeak
103       
104        Make/O/D/N=(numPeak,numTube) PeakTableX,peakTableY              //*2 to store x-location and peak height (y)
105        Edit peakTableX
106End
107
108Proc Identify_AllPeaks()
109
110        Variable ii,numTubes=8
111        String str="tube"
112       
113        ii=1
114        do
115                Identify_Peaks(str+num2str(ii),ii-1)
116                ii+=1
117        while(ii<=numTubes)
118
119End
120
121Proc Identify_Peaks(tubeStr,ind)
122        String tubeStr
123        Variable ind
124       
125        // must use a wave of pixels rather than "calculated" -- if calculated is used it only
126        // returns integer values for the peak locations
127       
128//      AutomaticallyFindPeaks() //-- this is a function that doesn't take any parameters - so
129// I need to pull this from the WM function to call the worker directly
130        Variable pBegin=0, pEnd= numpnts($(tubeStr))-1
131        Variable/C estimates= EstPeakNoiseAndSmfact($(tubeStr),pBegin, pEnd)
132        Variable noiselevel=real(estimates)
133        Variable smoothingFactor=imag(estimates)
134        Variable maxPeaks = 20
135        Variable minPeakPercent = 10
136       
137        AutoFindPeaksWorker($(tubeStr), $("tube_pixel"), pBegin, pEnd, maxPeaks, minPeakPercent, noiseLevel, smoothingFactor)
138// end WM function call
139
140        Sort WA_PeakCentersX WA_PeakCentersY,WA_PeakCentersX
141       
142        peakTableX[][ind] = WA_PeakCentersX[p]
143        peakTableY[][ind] = WA_PeakCentersY[p]
144       
145End
146
147// -- save a copy of the root:WA_PeakCentersY,root:WA_PeakCentersX values
148//    for later in case the fitting failed, then you can go back and re-do
149//
150// -- then plot:
151//
152//      Display peak_spacing_mm_ctr vs WA_PeakCentersX
153//
154//      Then do a "QuickFit" of the peak position to the data using a polynomial of order 3 (= quadratic)
155//
156// result is in W_coef, W_sigma
157//
158// -- an example of the "quickFit" command is below, so it can be programmed rather than the menu every time
159//¥Display peak_spacing_mm_ctr vs WA_PeakCentersX3
160//¥CurveFit/M=2/W=0/TBOX=(0x310) poly 3, peak_spacing_mm_ctr/X=WA_PeakCentersX3/D
161//  fit_peak_spacing_mm_ctr= poly(W_coef,x)
162//  W_coef={-571.42,1.1135,-4.2444e-05}
163//  V_chisq= 8.5841;V_npnts= 20;V_numNaNs= 0;V_numINFs= 0;
164//  V_startRow= 0;V_endRow= 19;
165//  W_sigma={0.595,0.00246,2.15e-06}
166//  Coefficient values ± one standard deviation
167//      K0      =-571.42 ± 0.595
168//      K1      =1.1135 ± 0.00246
169//      K2      =-4.2444e-05 ± 2.15e-06
170//
171//
172//
173// for (8) tubes, keep all of the fit coefficients
174//
175//¥make/O/D/N=(3,8) fit_coef
176//¥edit fit_coef
177//¥make/O/D/N=(3,8) fit_sigma
178//¥edit fit_sigma
179//
180// -- copy and paste in the W_coef and W_sigma values (or by a command)
181//
182
183
184Proc MakeTableForFitCoefs(numTube,numCoef)
185        Variable numTube,numCoef
186       
187        Make/O/D/N=(numCoef,numTube) TubeCoefTable,TubeSigmaTable               //
188        Edit TubeCoefTable
189End
190
191Proc PlotFit_AllPeaks()
192
193        Variable ii,numTubes=8
194       
195        ii=1
196        do
197                PlotFit_Peaks(ii-1)
198                ii+=1
199        while(ii<=numTubes)
200
201End
202
203Proc PlotFit_Peaks(ind)
204        Variable ind
205       
206        //hopefully 20 points - need better control of this
207        Duplicate/O WA_PeakCentersX, tmpX
208       
209        tmpX = peakTableX[p][ind]
210        Display peak_spacing_mm_ctr vs tmpX
211       
212        CurveFit/M=2/W=0/TBOX=(0x310) poly 3, peak_spacing_mm_ctr/X=tmpX/D
213       
214        TubeCoefTable[][ind] = W_coef[p]
215        TubeSigmaTable[][ind] = W_sigma[p]
216       
217End
218
219
220
221
222
223//¥Duplicate tube1 tube1_mm
224//¥tube1_mm = V_TubePix_to_mm(fit_coef[0][0],fit_coef[1][0],fit_coef[2][0],p)
225
226
227
228
229////////
230// then there are various display options:
231
232// adjust the center (pixel) of the tube:
233// - measCtr is the pixel location of the DEFINED "zero" peak
234// nominal Ctr is the pixel number of this DEFINED zer0 position, nominally nPix/2-1
235//
236// ( be sure to pick better names, and use a loop...)
237//      adj_tube = raw_tube[p+(measCtr-nominalCtr)]
238//
239// then fill and display a new matrix. The center will be reasonably well aligned, and will
240// get worse towards the ends of the tubes
241// (this may be the "preferred" way of displaying the raw data if the centers are far off)
242// -- this also may be what I need to start with to fit the data to locate the beam center.
243
244
245// I can also display the fully corrected tubes, where the y-axis is now in real-space mm
246// rather than arbitrary pixels. The x-axis is still tube nubmer.
247// -- do this with the procedure"
248//   Append_Adjusted_mm()               // name may change...
249//
250// -- the point of the appending is that it allows each tube to be plotted on an image plot
251// with its own y-axis. Every one will be different and will be non-linear. These conditions
252// BOTH prevent using any of the "normal" image plotting or manipulation routines.
253// - the gist is below:
254//
255//      duplicate tube1_mm adjusted_mm_edge
256//      InsertPoints 0,1, adjusted_mm_edge
257//      // be sure to use the correct set of coefficients
258//      adjusted_mm_edge[0] = V_TubePix_to_mm(fit_coef[0][0],fit_coef[1][0],fit_coef[2][0],-1)
259//     
260//      Display;AppendImage adjusted_pack vs {*,adjusted_mm_edge}
261
262
263Proc Make_mm_tubes()
264
265        Variable ii,numTubes=8
266       
267        ii=1
268        do
269                Duplicate $("tube"+num2str(ii)) $("tube"+num2str(ii)+"_mm")
270                $("tube"+num2str(ii)+"_mm") = V_TubePix_to_mm(TubeCoefTable[0][ii-1],TubeCoefTable[1][ii-1],TubeCoefTable[2][ii-1],p)
271                ii+=1
272        while(ii<=numTubes)
273       
274End
275
276
277Proc Append_Adjusted_mm()
278
279// a blank base image
280        Duplicate/O pack junk
281        junk=0
282        SetScale/I y -600,600,"", junk          // -600,600 is tooo large and not general
283        Display;Appendimage junk
284
285        Variable ii
286       
287        ii=1
288        do
289                make/O/D/N=(1,1127) $("tube"+num2str(ii)+"_mm_mat")=0   
290       
291                $("tube"+num2str(ii)+"_mm_mat")[0][] = $("tube"+num2str(ii))
292                SetScale/I x (ii-1),ii,"", $("tube"+num2str(ii)+"_mm_mat")              //builds up the x-axis
293               
294                duplicate/O $("tube"+num2str(ii)+"_mm") $("edge"+num2str(ii)+"_mm")
295                InsertPoints 0,1, $("edge"+num2str(ii)+"_mm")           //needs to be one point longer
296        // be sure to use the correct set of coefficients
297                $("edge"+num2str(ii)+"_mm")[0] = V_TubePix_to_mm(TubeCoefTable[0][0],TubeCoefTable[1][0],TubeCoefTable[2][0],-1)
298       
299                AppendImage $("tube"+num2str(ii)+"_mm_mat") vs {*,$("edge"+num2str(ii)+"_mm")}
300                ModifyImage $("tube"+num2str(ii)+"_mm_mat") ctab= {*,*,ColdWarm,0}
301       
302                ii+=1
303        while(ii < 9)
304       
305end
306
307
308
309////////////////////////////////
310
311Proc MakeMatrix_PixCenter()
312        Duplicate/O pack pack_centered
313       
314        Variable ii,numTubes=8
315       
316        ii=1
317        do
318                FillMatrix_Pix_Center(ii)
319                ii+=1
320        while(ii<=numTubes)
321       
322        Display;AppendImage pack_centered
323        ModifyImage pack_centered ctab= {*,*,ColdWarm,0}
324
325end
326
327//
328// this fills a matrix with the tubes center aligned ONLY, with the y-axis in pixels
329//
330// adj_tube = raw_tube[p+(measCtr-nominalCtr)]
331// finds the center automatically from the tubeN_mm wave, where it crosses zero
332//
333// Tube #1 is the "base" ans others are shifted to match that one's "zero"
334//
335// FindRoots/P=W_coef           can also be used to find the roots (but which one?)
336//
337Function FillMatrix_Pix_Center(ind)
338        Variable ind
339       
340        Wave adj=root:pack_centered
341        Wave tube1_mm = $("root:tube1_mm")
342        Wave tube = $("root:tube"+num2str(ind))
343        wave tube_mm = $("root:tube"+num2str(ind)+"_mm")
344
345        Variable base,shift,ii,num,pt
346       
347        num=numpnts(tube)
348       
349        FindLevel tube1_mm 0
350        base = round(V_LevelX)
351       
352       
353        FindLevel tube_mm 0
354        shift = round(V_LevelX)
355       
356        for(ii=0;ii<num;ii+=1)
357                pt = ii + (shift-base)
358                if(pt >= 0 && pt < num)
359                        adj[ind-1][ii] = tube[pt]
360                endif
361        endfor
362       
363        return(0)
364End
365
366
367// this fills a matrix with the tubes center aligned ONLY, with the y-axis in mm
368// -- there seems to be little reason to do this --
369// -- either keep pixels and align centers
370// -- OR -- use mm and append each tube with its own y-axis
371//
372Function FillAdjusted(ind)
373        Variable ind
374       
375        Wave adj=root:adjusted_pack
376        Wave tube1_mm
377        Wave tube = $("root:tube"+num2str(ind))
378        wave tube_mm = $("root:tube"+num2str(ind)+"_mm")
379
380        Variable base,shift,ii,num,pt
381       
382        num=numpnts(tube1_mm)
383       
384        FindLevel tube1_mm 0
385        base = round(V_LevelX)
386       
387       
388        FindLevel tube_mm 0
389        shift = round(V_LevelX)
390       
391        for(ii=0;ii<num;ii+=1)
392                pt = ii + (shift-base)
393                if(pt >= 0 && pt < num)
394                        adj[ind-1][ii] = tube[pt]
395                endif
396        endfor
397       
398        return(0)
399End
400
401
402// c0,c1,c2,pix
403// c0-c2 are the fit coefficients
404// pix is the test pixel
405//
406// returns the distance in mm (relative to ctr pixel)
407// ctr is the center pixel, as defined when fitting to quadratic
408//
409Function V_TubePix_to_mm(c0,c1,c2,pix)
410        Variable c0,c1,c2,pix
411       
412        Variable dist
413        dist = c0 + c1*pix + c2*pix*pix
414       
415        return(dist)
416End
417
418////
419
420
421
422// set the (linear) range of the y-axis of the matrix to be the
423// range of the 1st tube. This is completely arbitrary
424//
425Proc Interpolate_mm_tubes()
426
427        Duplicate/O pack pack_image
428
429        Variable ii,numTubes=8
430        Variable p1,p2
431        p1 = tube1_mm[0]
432        p2 = tube1_mm[numpnts(tube1_mm)-1]
433       
434        SetScale/I y p1,p2,"", pack_image
435       
436        // then make a temporary 1D wave to help with the interpolation
437        Duplicate/O tube1_mm lin_mm,lin_val
438        SetScale/I x p1,p2,"", lin_mm
439        lin_mm = x                      //fill with the linear mm spacing
440        lin_val=0
441       
442        ii=1
443        do
444                lin_val = interp(lin_mm, $("tube"+num2str(ii)+"_mm"), $("tube"+num2str(ii)))
445                pack_image[ii-1][] = lin_val[q]
446               
447                ii+=1
448        while(ii<=numTubes)
449       
450        display;appendimage pack_image
451        ModifyImage pack_image ctab= {*,*,ColdWarm,0}
452       
453End
454
455
456
457
458
459
460
461
462
463
464
465
466
467// this doesn't work - the interpolation step doesn't do what I want
468// and the plot of the triplet with f(z) for color doesn't fill space like I want
469Proc AnotherExample()
470
471        Concatenate/O/NP {tube1_mm,tube2_mm,tube3_mm,tube4_mm,tube5_mm,tube6_mm,tube7_mm,tube8_mm},cat_mm
472        Concatenate/O/NP {tube1,tube2,tube3,tube4,tube5,tube6,tube7,tube8},cat_tubes
473        Duplicate/O cat_mm,cat_num
474        Variable num=1127
475        cat_num[0,num-1]=1
476        cat_num[num,2*num-1]=2
477        cat_num[2*num,3*num-1]=3
478        cat_num[3*num,4*num-1]=4
479        cat_num[4*num,5*num-1]=5
480        cat_num[5*num,6*num-1]=6
481        cat_num[6*num,7*num-1]=7
482        cat_num[7*num,8*num-1]=8
483
484        Display cat_mm vs cat_num
485        ModifyGraph mode=3,marker=9
486        ModifyGraph zColor(cat_mm)={cat_tubes,*,*,ColdWarm,0}
487
488        Concatenate/O {cat_num,cat_mm,cat_tubes}, tripletWave
489        ImageInterpolate Kriging tripletWave
490        AppendImage M_InterpolatedImage
491
492//      Make/O/N=20 xWave=enoise(4),yWave=enoise(5),zWave=enoise(6)  // Random points
493//      Display yWave vs xWave
494//      ModifyGraph mode=3,marker=19
495//      ModifyGraph zColor(yWave)={zWave,*,*,Rainbow,0}
496//
497//      Concatenate/O {xWave,yWave,zWave}, tripletWave
498//      ImageInterpolate/S={-5,0.1,5,-5,0.1,5} voronoi tripletWave
499//      AppendImage M_InterpolatedImage
500
501end
502
503// this desn't work either...
504// (same y-axis for the entire image, which is not the case for the tubes)
505//
506// from the WM help file:
507// Plotting a 2D Z Wave With 1D X and Y Center Data
508//
509Function MakeEdgesWave(centers, edgesWave)
510        Wave centers    // Input
511        Wave edgesWave  // Receives output
512       
513        Variable N=numpnts(centers)
514        Redimension/N=(N+1) edgesWave
515
516        edgesWave[0]=centers[0]-0.5*(centers[1]-centers[0])
517        edgesWave[N]=centers[N-1]+0.5*(centers[N-1]-centers[N-2])
518        edgesWave[1,N-1]=centers[p]-0.5*(centers[p]-centers[p-1])
519End
520
521//This function demonstrates the use of MakeEdgesWave:
522Function DemoPlotXYZAsImage()
523        Make/O mat={{0,1,2},{2,3,4},{3,4,5}}    // Matrix containing Z values
524        Make/O centersX = {1, 2.5, 5}           // X centers wave
525        Make/O centersY = {300, 400, 600}               // Y centers wave
526        Make/O edgesX; MakeEdgesWave(centersX, edgesX)  // Create X edges wave
527        Make/O edgesY; MakeEdgesWave(centersY, edgesY)  // Create Y edges wave
528        Display; AppendImage mat vs {edgesX,edgesY}
529End
530
Note: See TracBrowser for help on using the repository browser.