Ignore:
Timestamp:
Dec 2, 2013 1:36:53 PM (9 years ago)
Author:
srkline
Message:

Added utilities to RealSpaceModeling? to allow previously saved matrices to be added together. currently they are strictly added - that is overlapping voxels are added together, which may not be the optimal behavior, but this can be changed as needed.

Also added utility to allow arbitrary rotation of a 3D object. No interpolation is done, so the rotation is lossy. Rotation is applied in xyz order, which can make a difference. Beware.

Documentation to follow (eventually).

Location:
sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Cubes.ipf

    r925 r926  
    620620         
    621621        Qmstr = num2str(FFT_QMaxReal-2*deltaQ)          //using FFT_QMaxReal gives NaN at edges of interp2Dslice 
    622         str = "FakeQxQy(0,"+qmstr+",-"+qmstr+","+qmstr+","+num2str(num)+",\"fakeHalf\",1,0)"            //overwrites the data folder, no plot 
     622        str = "FakeQxQy(0,"+qmstr+",-"+qmstr+","+qmstr+","+num2str(num)+",\"USANS_Half\",1,0)"          //overwrites the data folder, no plot 
    623623        Execute str 
    624624 
     
    628628        FFT_Get2DSlice("") 
    629629// interpolate to the USANS-sized detector 
    630         Interpolate2DSliceToData("fakeHalf") 
     630        Interpolate2DSliceToData("USANS_Half") 
    631631        Wave interp2Dslice = $("root:interp2Dslice") 
    632632 
     
    642642 
    643643        FFT_aUSANS_i *= 4               //why the factor of 4??? 
     644         
     645 
     646        Execute "USANSDetectorHalf()" 
     647 
    644648         
    645649        // not needed -- occupancy and the normal FFT scaling are already done -- the slice is already on absolute scale 
     
    677681Function Isotropic_FFT_to_USANS() 
    678682 
    679 ////////////// 
    680683        Variable maxQx,Qy,num,deltaQ 
    681684        String str,qmStr 
     
    696699         
    697700        Qmstr = num2str(FFT_QMaxReal-2*deltaQ)          //using FFT_QMaxReal gives NaN at edges of interp2Dslice 
    698         str = "FakeQxQy(0,"+qmstr+",-"+qmstr+","+qmstr+","+num2str(num)+",\"fakeHalf\",1,0)"            //overwrites the data folder, no plot 
     701        str = "FakeQxQy(0,"+qmstr+",-"+qmstr+","+qmstr+","+num2str(num)+",\"USANS_Half\",1,0)"          //overwrites the data folder, no plot 
    699702        Execute str 
    700703 
     
    734737         
    735738// interpolate to the USANS-sized detector, this is equivalent to: 
    736 //      Interpolate2DSliceToData("fakeHalf")            //except using avg2d rather than "calc = detPlane" (the anisotropic result) 
    737         String folderStr = "fakeHalf" 
     739//      Interpolate2DSliceToData("USANS_Half")          //except using avg2d rather than "calc = detPlane" (the anisotropic result) 
     740        String folderStr = "USANS_Half" 
    738741        Wave data = $("root:"+folderStr+":"+folderStr+"_lin") 
    739742         
     
    766769        Execute "DisplayFFT_to2DAverage()" 
    767770 
     771        Execute "USANSDetectorHalf()" 
     772 
    768773// if you want to scale this to match a model calculation, you'll need the volume fraction 
    769774        Variable occ = VolumeFraction_Occ(root:mat) 
     
    771776                 
    772777        return(0) 
     778End 
     779 
     780Proc USANSDetectorHalf() 
     781        DoWindow USANS_DetectorHalf 
     782        if(V_flag == 0) 
     783                PauseUpdate; Silent 1           // building window... 
     784                Display /W=(898,808,1083,1150) 
     785                DoWindow/C USANS_DetectorHalf 
     786                AppendImage/T interp2DSlice_log 
     787                ModifyImage interp2DSlice_log ctab= {*,*,YellowHot,0} 
     788                ModifyGraph margin(left)=14,margin(bottom)=14,margin(top)=14,margin(right)=14,height={Aspect,2} 
     789                ModifyGraph mirror=2 
     790                ModifyGraph nticks=4 
     791                ModifyGraph minor=1 
     792                ModifyGraph fSize=9 
     793                ModifyGraph standoff=0 
     794                ModifyGraph tkLblRot(left)=90 
     795                ModifyGraph btLen=3 
     796                ModifyGraph tlOffset=-2 
     797                SetAxis/A/R left 
     798        endif 
    773799End 
    774800 
  • sans/Dev/trunk/NCNR_User_Procedures/Analysis/Alpha/Tinker/FFT_Panel.ipf

    r925 r926  
    9191        Button FFTButton_5,pos={13,218},size={120,20},proc=FFTDrawZCylinderButtonProc,title="Draw Cylinder" 
    9292//      Button FFTButton_6,pos={134,79},size={90,20},proc=FFTEraseMatrixButtonProc,title="Erase Matrix" 
    93         Button FFTButton_6a,pos={160,79},size={60,20},proc=FFTSaveMatrixButtonProc,title="Save" 
    94         Button FFTButton_6b,pos={240,79},size={60,20},proc=FFTLoadMatrixButtonProc,title="Load" 
     93        Button FFTButton_6a,pos={140,79},size={50,20},proc=FFTSaveMatrixButtonProc,title="Save" 
     94        Button FFTButton_6b,pos={200,79},size={50,20},proc=FFTLoadMatrixButtonProc,title="Load" 
     95        Button FFTButton_6c,pos={260,79},size={50,20},proc=FFT_AddMatrixButtonProc,title="Add" 
    9596        Button FFTButton_7,pos={13,329},size={130,20},proc=FFT_BinnedSpheresButtonProc,title="Do Binned Debye" 
    9697        Button FFTButton_7a,pos={180,329},size={130,20},proc=FFT_PlotResultsButtonProc,title="Plot Binned Results" 
     
    147148        sprintf str,"FFT_T=%g;FFT_N=%d;FFT_SolventSLD=%d;",FFT_T,FFT_N,FFT_SolventSLD 
    148149        Note mat,str 
    149         Save/C/P=home mat as fileStr    //will ask for a file name, save as Igor Binary 
     150        Save/C/I/P=home mat as fileStr  //will ask for a file name, save as Igor Binary 
    150151        Note/K mat                      //kill wave note on exiting since I don't properly update this anywhere else 
    151152                         
     
    173174End 
    174175 
    175  
     176// /H flag on the LoadWave command severs the connection with the binary file 
     177// -- this seems to be important - otherwise I get odd results and the wave (on disk) can change! 
     178// 
    176179Function ReloadMatrix(fileStr) 
    177180        String fileStr 
    178181         
    179                 LoadWave/M/O/W/P=home           fileStr         //will ask for a file, Igor Binary format is assumed here 
    180                 if(V_flag == 0) 
    181                         return(0)               //user cancel 
     182        LoadWave/H/M/O/W/P=home         fileStr         //will ask for a file, Igor Binary format is assumed here 
     183        if(V_flag == 0) 
     184                return(0)               //user cancel 
     185        endif 
     186         
     187        String str 
     188        str=note(mat) 
     189        NVAR FFT_T = root:FFT_T 
     190        NVAR FFT_N = root:FFT_N 
     191        NVAR FFT_SolventSLD = root:FFT_SolventSLD 
     192         
     193        FFT_T = NumberByKey("FFT_T", str, "=" ,";") 
     194        FFT_N = NumberByKey("FFT_N", str, "=" ,";") 
     195        FFT_SolventSLD = NumberByKey("FFT_SolventSLD", str, "=" ,";") 
     196 
     197// if I got bad values, put in default values                    
     198        if(numtype(FFT_T) != 0 ) 
     199                FFT_T = 5 
     200        endif 
     201        if(numtype(FFT_N) != 0 ) 
     202                FFT_N = DimSize(mat,0) 
     203        endif 
     204        if(numtype(FFT_SolventSLD) != 0 ) 
     205                FFT_SolventSLD = 0 
     206        endif                    
     207 
     208        ColorizeGizmo() 
     209         
     210        Print "Loaded matrix parameters = ",str 
     211        Execute "NumberOfPoints()" 
     212                         
     213        return(0) 
     214end 
     215 
     216// load in a previously saved matrix, and reset FFT_N, FFT_T and solvent 
     217// from the wave note when saved 
     218Function FFT_AddMatrixButtonProc(ba) : ButtonControl 
     219        STRUCT WMButtonAction &ba 
     220         
     221        String win = ba.win 
     222 
     223        switch (ba.eventCode) 
     224                case 2: 
     225                        // click code here 
     226                        String fileStr="" 
     227                        LoadAndAddMatrix(fileStr) 
     228                         
     229                        break 
     230        endswitch 
     231 
     232        return 0 
     233End 
     234 
     235// This function will load and add the matrix to whatever is currently present 
     236// - it is checked that the N, T, and SolventSLD are the same 
     237// 
     238// - it currently SUMS the voxels - so if there is overlap, you get a new value 
     239// 
     240Function LoadAndAddMatrix(fileStr) 
     241        String fileStr 
     242 
     243        String toLoadStr="" 
     244 
     245// if the passed in file name is null, pick a file       
     246        if(strlen(fileStr)==0) 
     247                toLoadStr = DoOpenFileDialog("Pick the binary file to load") 
     248                if(strlen(toLoadStr)==0) 
     249                        return(0) 
    182250                endif 
     251        endif 
     252        fileStr=toLoadStr 
     253 
     254        // make sure that N and T are correct, just read in the note 
     255        String noteStr=LoadNoteFunc("home",fileStr) 
     256        String abortStr 
     257         
     258        // current values 
     259        NVAR FFT_T = root:FFT_T 
     260        NVAR FFT_N = root:FFT_N 
     261        NVAR FFT_SolventSLD = root:FFT_SolventSLD 
     262 
     263// possible set to load 
     264        Variable test_T,test_N,test_SolventSLD 
     265        test_T = NumberByKey("FFT_T", noteStr, "=" ,";") 
     266        test_N = NumberByKey("FFT_N", noteStr, "=" ,";") 
     267        test_SolventSLD = NumberByKey("FFT_SolventSLD", noteStr, "=" ,";") 
     268         
     269// if I got bad values, warn and abort                   
     270        if(FFT_T != test_T) 
     271                abortStr = "Current T = "+num2str(FFT_T)+" and Selected T = "+num2str(test_T)+", aborting load" 
     272                Abort abortStr 
     273        endif 
     274        if(FFT_N != test_N) 
     275                abortStr = "Current N = "+num2str(FFT_N)+" and Selected N = "+num2str(test_N)+", aborting load" 
     276                Abort abortStr 
     277        endif 
     278        if(FFT_SolventSLD != test_SolventSLD) 
     279                abortStr = "Current SolventSLD = "+num2str(FFT_SolventSLD)+" and Selected SolventSLD = "+num2str(test_SolventSLD)+", aborting load" 
     280                Abort abortStr 
     281        endif    
    183282                 
    184                 String str 
    185                 str=note(mat) 
    186                 NVAR FFT_T = root:FFT_T 
    187                 NVAR FFT_N = root:FFT_N 
    188                 NVAR FFT_SolventSLD = root:FFT_SolventSLD 
     283        // OK, then load and add the matrix 
     284        // the loaded matrix is "mat" as usual, so make a copy of the current first 
     285        Duplicate/O mat tmpMat 
     286        LoadWave/H/M/O/W/P=home         fileStr         //will ask for a file, Igor Binary format is assumed here 
     287        if(V_flag == 0) 
     288                return(0)               //user cancel 
     289        endif 
     290         
     291        Wave mat=root:mat 
     292        FastOp mat = mat + tmpMat 
     293         
     294        KillWaves/Z tmpMat 
     295         
     296        Print "Loaded matrix parameters = ",noteStr 
     297        Execute "NumberOfPoints()" 
     298         
     299        ColorizeGizmo() 
    189300                 
    190                 FFT_T = NumberByKey("FFT_T", str, "=" ,";") 
    191                 FFT_N = NumberByKey("FFT_N", str, "=" ,";") 
    192                 FFT_SolventSLD = NumberByKey("FFT_SolventSLD", str, "=" ,";") 
    193  
    194 // if I got bad values, put in default values                    
    195                 if(numtype(FFT_T) != 0 ) 
    196                         FFT_T = 5 
    197                 endif 
    198                 if(numtype(FFT_N) != 0 ) 
    199                         FFT_N = DimSize(mat,0) 
    200                 endif 
    201                 if(numtype(FFT_SolventSLD) != 0 ) 
    202                         FFT_SolventSLD = 0 
    203                 endif                    
    204                  
    205                 Print "Loaded matrix parameters = ",str 
    206                 Execute "NumberOfPoints()" 
    207                          
    208301        return(0) 
    209302end 
     
    446539        String ctrlName 
    447540         
    448         Variable degree=45,sense=1 
    449         Prompt degree, "Degrees of rotation around Z-axis:" 
    450 //      Prompt sense, "Direction of rotation:",popup,"CW;CCW;" 
    451         DoPrompt "Enter parameters for rotation", degree 
     541        Variable angleX=45,angleY=0,angleZ=0 
     542        Prompt angleX, "Degrees of rotation around X-axis:" 
     543        Prompt angleY, "Degrees of rotation around Y-axis:" 
     544        Prompt angleZ, "Degrees of rotation around Z-axis:" 
     545        DoPrompt "Enter angles for rotation", angleX,angleY,angleZ 
    452546         
    453547        if (V_Flag) 
    454548                return 0                                                                        // user canceled 
    455549        endif 
    456          
    457         fFFT_RotateMat(degree) 
     550 
     551        XYZRotate(angleX,angleY,angleZ) 
     552 
     553 
     554//////////// old way, only one angle 
     555//      Variable degree=45,sense=1 
     556//      Prompt degree, "Degrees of rotation around Z-axis:" 
     557////    Prompt sense, "Direction of rotation:",popup,"CW;CCW;" 
     558//      DoPrompt "Enter parameters for rotation", degree 
     559//       
     560//      if (V_Flag) 
     561//              return 0                                                                        // user canceled 
     562//      endif 
     563         
     564 
     565        // old way using ImageRotate that interpolates, and is only around Z-axis 
     566//      fFFT_RotateMat(degree) 
    458567        return(0) 
    459568End 
     
    692801        Wave mat=root:mat 
    693802        NVAR val=root:FFT_SolventSLD 
    694         mat=val 
     803        FastOp mat=(val) 
    695804        return(0) 
    696805End 
     
    779888         
    780889End 
    781  
    782  
    783  
    784890 
    785891 
     
    9251031End 
    9261032 
    927  
     1033////////////// my functions to rotate the matrix in a XYZ coordinate system 
     1034// 
     1035// 
     1036// definitely not generic, is expecting NxNxN volume 
     1037// 
     1038// not very friendly in that it "clips" anything that rotates out of the volume 
     1039// 
     1040// does translate the center of the box to 000, rotates, then translates back 
     1041// 
     1042// friendly in the sense that the rotated matrix is the same size as the original. 
     1043//  --this is important for my final application (FFT) 
     1044// 
     1045// does no interpolation of values, so be sure to keep a copy of the original 
     1046// -- multiple rotation steps are going to make a mess of things. 
     1047// 
     1048// The multi axis rotation is done as one step, and probably violates every conventional 
     1049//  coordinate system. The rotation is applied as RxRyRz, but this could easily be changed 
     1050// 
     1051// I just want it to be correct, so speed was not an issue. 
     1052// -- it's nested for loops. 
     1053// -- it's working with the full matrix, even when 99% is empty. 
     1054// 
     1055// 20 NOV 2013 SRK 
     1056// 
     1057// 
     1058 
     1059// 
     1060// mat is the input volume 
     1061// rotVol is the output rotated volume 
     1062Function XYZRotate(angleX,angleY,angleZ) 
     1063        Variable angleX,angleY,angleZ 
     1064         
     1065        NVAR FFT_N=root:FFT_N 
     1066        WAVE mat=root:mat 
     1067        Variable dist=FFT_N/2 
     1068 
     1069 
     1070// convert the NxNxN into 3xN xyz locations + wave of "w" values named "values" 
     1071        fVolumeToXYZTriplet(mat,"trip") 
     1072        Wave trip=root:trip 
     1073 
     1074// translate to get the center of the xyz values to 0,0,0        
     1075        fTranslateCoordinate(trip,dist)         //subtracts dist 
     1076 
     1077// do the rotation as a matrix multiplication    
     1078// putting zero is no rotation around that axis 
     1079        DoRotation(trip,angleX,angleY,angleZ) 
     1080        Wave rotated=root:rotated 
     1081 
     1082// translate back to a 0->N based coordinate 
     1083        fTranslateCoordinate(rotated,-dist) 
     1084        Wave values=root:values 
     1085 
     1086// convert the triplet back to a volume 
     1087// this CLIPS anything that has rotated out of the NxNxN volume 
     1088        fXYZTripletToVolume(rotated,values,"rotVol",FFT_N) 
     1089 
     1090// clean up by killng the extra waves that were generated 
     1091// 
     1092        KillWaves/Z trip,rotated,values 
     1093         
     1094        Wave rotVol=root:rotVol 
     1095        mat=rotVol 
     1096         
     1097        return(0) 
     1098End 
     1099 
     1100 
     1101 
     1102Function fVolumeToXYZTriplet(matrixWave, outputName) 
     1103        Wave matrixWave 
     1104        String outputName        
     1105         
     1106        Variable dimx=DimSize(matrixWave,0) 
     1107        Variable dimy=DimSize(matrixWave,1) 
     1108        Variable dimz=DimSize(matrixWave,2) 
     1109        Variable rows=dimx*dimy*dimz 
     1110        Make/O/N=(3,rows) $outputName 
     1111        Make/O/N=(rows) values 
     1112        WAVE TripletWave= $outputName 
     1113        Wave values=values 
     1114         
     1115 
     1116        Variable ii,jj,kk,count=0 
     1117        Variable xVal,yVal,zval 
     1118        for(kk=0;kk<dimz;kk+=1)                 // kk is z (layer) 
     1119                zval=kk 
     1120                for(jj=0;jj<dimy;jj+=1)         // jj is y (column) 
     1121                        yVal=jj 
     1122                        for(ii=0;ii<dimx;ii+=1) // ii is x (row) 
     1123                                xVal=ii 
     1124                                TripletWave[0][count]=xVal 
     1125                                TripletWave[1][count]=yVal 
     1126                                TripletWave[2][count]=zval 
     1127                                values[count]=matrixWave[ii][jj][kk] // value at [row][col][lay] 
     1128                                count+=1 
     1129                        endfor 
     1130                endfor 
     1131        endfor 
     1132         
     1133        return(0) 
     1134End 
     1135 
     1136 
     1137Function fXYZTripletToVolume(triplet, values, outputName, outputDim) 
     1138        Wave triplet,values 
     1139        String outputName 
     1140        Variable outputDim 
     1141         
     1142        Variable numPt=DimSize(triplet,1) 
     1143 
     1144        Variable num = outputDim 
     1145         
     1146        Make/O/B/N=(num,num,num) $outputName 
     1147        WAVE newVol= $outputName 
     1148 
     1149        FastOp newVol = 0 
     1150         
     1151        Variable ii,jj,kk,count=0 
     1152        Variable xVal,yVal,zval 
     1153        Variable xOK, yOK, zOK 
     1154 
     1155         
     1156        for(ii=0;ii<numPt;ii+=1) 
     1157                xval = round(triplet[0][ii]) 
     1158                yval = round(triplet[1][ii]) 
     1159                zval = round(triplet[2][ii]) 
     1160                 
     1161                // round and keep in bounds (returns truth) 
     1162                xOK = inRange(xval,0,num-1) 
     1163                yOK = inRange(yval,0,num-1) 
     1164                zOK = inRange(zval,0,num-1) 
     1165                 
     1166                if(xOK && yOK && zOK) 
     1167                        newVol[xval][yval][zval] = values[ii] 
     1168                endif 
     1169                         
     1170        endfor 
     1171 
     1172 
     1173        return(0) 
     1174End 
     1175 
     1176// if val < lo or > hi, bad val 
     1177// fi both of these pass, pt is OK 
     1178ThreadSafe Static Function inRange(val, lo, hi) 
     1179        Variable val, lo, hi 
     1180  
     1181        if(val < lo) 
     1182                return (0) 
     1183        endif 
     1184        if(val > hi) 
     1185                return (0) 
     1186        endif 
     1187         
     1188        return(1) 
     1189  
     1190End 
     1191 
     1192// Rotation is applied in the order Rx Ry Rz 
     1193Function DoRotation(triplet,angleX,angleY,angleZ) 
     1194        Wave triplet 
     1195        Variable angleX,angleY,angleZ 
     1196         
     1197        Variable thetaX,thetaY,thetaZ 
     1198        thetaX = angleX*pi/180          // convert degrees to radians 
     1199        thetaY = angleY*pi/180          // convert degrees to radians 
     1200        thetaZ = angleZ*pi/180          // convert degrees to radians 
     1201         
     1202        Make/O/D/N=(3,3) Rx,Ry,Rz 
     1203        Rx=0 
     1204        Ry=0 
     1205        Rz=0 
     1206         
     1207        Rx[0][0] = 1 
     1208        Rx[1][1] = cos(thetaX) 
     1209        Rx[1][2] = -sin(thetaX) 
     1210        Rx[2][1] = sin(thetaX) 
     1211        Rx[2][2] = cos(thetaX) 
     1212         
     1213        Ry[0][0] = cos(thetaY) 
     1214        Ry[0][2] = sin(thetaY) 
     1215        Ry[1][1] = 1 
     1216        Ry[2][0] = -sin(thetaY) 
     1217        Ry[2][2] = cos(thetaY) 
     1218         
     1219        Rz[0][0] = cos(thetaZ) 
     1220        Rz[0][1] = -sin(thetaZ) 
     1221        Rz[1][0] = sin(thetaZ) 
     1222        Rz[1][1] = cos(thetaZ) 
     1223        Rz[2][2] = 1     
     1224         
     1225         
     1226        MatrixOp/O rotated = Rx x Ry x Rz x triplet 
     1227         
     1228         
     1229        return(0) 
     1230end 
     1231 
     1232 
     1233// the rotation matrix, as I copied from wikipedia, is a rotation around (0,0,0), not 
     1234// the center of the gizmo plot 
     1235Function fTranslateCoordinate(trip,dist) 
     1236        Wave trip 
     1237        Variable dist 
     1238         
     1239        MatrixOp/O trip = trip - dist 
     1240         
     1241        return(0) 
     1242End 
     1243 
Note: See TracChangeset for help on using the changeset viewer.